Unraveling the spatiotemporal brain dynamics during a simulated reach-to-eat task

The reach-to-eat task involves a sequence of action components including looking, reaching, grasping, and feeding. While cortical representations of individual action components have been mapped in human functional magnetic resonance imaging (fMRI) studies, little is known about the continuous spatiotemporal dynamics among these representations during the reach-to-eat task. In a periodic event-related fMRI experiment, subjects were scanned while they reached toward a food image, grasped the virtual food, and brought it to their mouth within each 16-s cycle. Fourier-based analysis of fMRI time series revealed periodic signals and noise distributed across the brain. Independent component analysis was used to remove periodic or aperiodic motion artifacts. Time-frequency analysis was used to analyze the temporal characteristics of periodic signals in each voxel. Circular statistics was then used to estimate mean phase angles of periodic signals and select voxels based on the distribution of phase angles. By sorting mean phase angles across regions, we were able to show the real-time spatiotemporal brain dynamics as continuous traveling waves over the cortical surface. The activation sequence consisted of approximately the following stages: (1) stimulus related activations in occipital and temporal cortices; (2) movement planning related activations in dorsal premotor and superior parietal cortices; (3) reaching related activations in primary sensorimotor cortex and supplementary motor area; (4) grasping related activations in postcentral gyrus and sulcus; (5) feeding related activations in orofacial areas. These results suggest that phase-encoded design and analysis can be used to unravel sequential activations among brain regions during a simulated reach-to-eat task.


Introduction
Eating with the hand and arm is one of the many daily activities that take place in the peripersonal space. Even a trivial task such as eating popcorn involves precise coordination among sight, touch, movement, and proprioception. A typical reach-to-eat task consists of a sequence of action components: (1) looking at the food; (2) reaching toward the food; (3) grasping the food; and (4) transporting the food to the mouth (Flindall and Gonzalez, 2013;Foroud and Whishaw, 2012;Quinlan and Culham, 2015;Sacrey et al., 2012;Sveistrup et al., 2008). Over the last two decades, human functional magnetic resonance imaging (fMRI) studies have revealed frontal and parietal regions that are involved in visually-guided actions, such as saccades, pointing, reaching, and grasping (see review in Culham and Valyear, 2006;Filimon, 2010;Gallivan and Culham, 2015;Grafton, 2010;Huang and Sereno, 2018;Sereno and Huang, 2014;Tunik et al., 2007;Vesia and Crawford, 2012). Furthermore, sensorimotor representations of hands (fingers), face, lips, jaw, and tongue were mapped in fMRI experiments using passive tactile stimulation or voluntary movements (Carey et al., 2017;Chen et al., 2017;Grabski et al., 2012;Huang and Sereno, 2007;Huang et al., 2012Huang et al., , 2017Meier et al., 2008;Miyamoto et al., 2006). While cortical visual, tactile, motor, and multisensory representations supporting visually-guided actions have been mapped individually, little is known about the spatiotemporal dynamics among these representations during a continuous reach-to-eat task.
To map cortical representations of individual action components in fMRI experiments, previous studies have used randomized block designs where each block consisted of the same action repeated multiple times (e.g., Filimon et al., 2007Filimon et al., , 2009)or event-related designswhere a single action or an action sequence constituted an event (Cavina-Pratesi et al., 2010Gallivan et al., 2011Gallivan et al., , 2016Rossit et al., 2013). Functional images are typically analyzed using the general linear model (GLM) that computes the overall fit between the design matrix and the fMRI time series of each voxel (Ward, 2002). Brain regions specialized for a specific action component or sequence are then identified and differentiated from the others using a subtraction method (i.e., linear contrasts between conditions) and by comparing the magnitudes of peak hemodynamic response within a time window after the event onset. However, the resulting "static" maps only show distinct brain regions containing statistically significant activations without revealing the temporal sequence among them.
The methods just described do not take full advantage of the temporal resolution of standard fMRI, which at approximately 1 s, is good enough to begin to resolve differences in activation latency among brain regions underlying cognitive, sensory, and/or motor processes in naturalistic tasks (Ariani et al., 2018;Cunnington et al., 2005;Formisano and Goebel, 2003;Formisano et al., 2002;Menon et al., 1998;Richter et al., 1997;Smolders et al., 2007;Sun et al., 2005;Weilke et al., 2001;Wildgruber et al., 1997;Windischberger et al., 2008). To map the topological organization of primary and association cortices, the phase-encoded design has been developed to make full use of the temporal precision of fMRI signals (Besle et al., 2013;Engel, 2012;Engel et al., 1994;Harvey et al., 2013;Sereno, 2007, 2013;Overduin and Servos, 2004;Sereno and Huang, 2006;Sereno et al., 1995Sereno et al., , 2001Talavage et al., 2004). It has also been used to map somatotopic organization of motor cortex by performing a sequence of movements periodically (Carey et al., 2017;Sood and Sereno, 2016). In these studies, the Fourier transform was performed on the fMRI time series of each voxel. The resulting power spectrum was used to measure activation significance, and then the phase angle at the stimulus (or task) frequency was used to measure the temporal offset (delay) of significant periodic activations. The temporal sequence of activations among brain regions was then visualized by animating successive iso-phase contours as traveling waves over the cortical surface (Engel et al., 1994;Sereno and Huang, 2006;Sereno et al., 1995;Windischberger et al., 2008). The surface-based traveling wave approach facilitates understanding the hemodynamic response as both temporal and spatial processes (Aquino et al., 2012;Lacy et al., 2016).
In this study, we explored the use of phase-encoded design and analysis for mapping the spatiotemporal dynamics among brain regions during a simulated reach-to-eat task, where subjects periodically and sequentially reached toward a food image, grasped the virtual food, and brought it to their mouth within each 16-s cycle. The design of the simulated eating task allows experiments to be performed in the confines of an MRI scanner where safety and sanitary rules preclude the presence of real food while capturing the coordinated, sequential movements of an actual eating task. It is expected that periodic reach-to-eat movements will result in periodic activations in cortical representations of action components individually mapped by GLM-based methods in previous studies (e.g., Cavina-Pratesi et al., 2010Gallivan et al., 2011Gallivan et al., , 2016Rossit et al., 2013). Given that fMRI has limited temporal resolution, we aim to determine whether it is possible to further resolve the relative activation latencies (i.e., early, intermediate, and late phases) among brain regions containing periodic activations. Toward this end, we used time-frequency analysis to analyze the temporal characteristics of fMRI time series . Circular statistics was then used to estimate mean phase angles of periodic signals and select voxels based on the distribution of phase angles. Mean phase angles were sorted and visualized as continuous traveling waves over the cortical surface, which revealed complex spatiotemporal brain dynamics during a compound action consisting of looking, reaching, grasping, and feeding components. In addition to the directions of propagation, these traveling waves revealed the interaction and coherence among brain regions, which are not obtainable by mapping each action component in isolation.

Participants
Nine healthy right-hand dominant subjects (18-30 years; 3 males, 6 females) with normal or corrected-to-normal vision participated in this study. All subjects gave written informed consent according to protocols approved by the Human Research Protections Program of the University of California, San Diego (UCSD).

Experimental design and setup
Each subject participated in a brief training session in a custom-built mock MRI scanner before an actual fMRI session. In both sessions, the subject lay supine with head tilted forward inside a head coil, where the head was supported and constrained by foam padding. A semitranslucent screen was mounted~15 cm in front of the face, allowing the subject to directly view and touch the screen. Visual stimuli were programmed in the C language using the OpenGL Performer library (Huang et al., 2015), rendered at a resolution of 1024 Â 768 pixels on a LCD projector (Dell 3300 MP), and back-projected onto a region of 35 Â 26 cm (~98 Â 82 field of view; Fig. 1). Throughout each session, the subject fixated a central cross against a black background while attending to an image of snacks (size ¼~28 Â 28 ; center ¼~26 eccentricity) located in the lower right visual field. There were 16 periodic and evenly-spaced trials (cycles) in each 256-s run. At the beginning of each run, a food image appeared against a black background, which signaled the onset (0 s) of the first trial. The image stayed on the screen for 16 s and was then replaced by a different image randomly selected from six ( Fig. 1), which signaled the onset of the next trial. Immediately after the detection of an image swap, the subject slowly reached toward the food image with the right forearm, touched and grasped the virtual food on the screen (with the thumb opposing the index, middle, and ring fingers), brought it into the mouth by touching both the upper and lower lips with finger tips, and then rested the hand under the chin. A sequence of reach-to-eat movements was performed continuously at a self-paced speed for 10-12 s, followed by a brief rest, within each 16-s trial (cycle). The subject was informed that the task was similar to the scenario Fig. 1. Stimulus layout and food images presented during the reach-to-eat task. At the onset of each 16-s cycle, a different image was randomly selected from six and shown in the lower right visual field for the entire cycle. of slowly eating popcorn from a container near the lower face while watching a movie straight ahead (without directly looking at the popcorn).
Head movement related artifacts in functional images are a major concern for most sensorimotor fMRI experiments. In previous studies, a number of preventive measures have been attempted to minimize head movements. First, subjects can be trained to remain completely still while performing motor tasks within a mock MRI scanner before they participate in actual fMRI experiments. Second, head movements can be reduced by bracing the shoulder and upper arm while allowing restricted forearm movements away from the head (Cavina-Pratesi et al., 2010Gallivan et al., 2011Gallivan et al., , 2016Rossit et al., 2013). Third, the head can be immobilized using foam pads, thermoplastic sheets, and/or a bite bar made with individual dental impression Filimon et al., 2007Filimon et al., , 2009Huang et al., 2012Huang et al., , 2017. In the present study, we recruited experienced subjects who had previously participated in multiple fMRI experiments and had been trained to stay motionless in mock and real MRI scanners. To clear the pathway for hand-to-mouth movements, we did not use thermoplastic sheets or bite bars. Instead, we used tight foam padding to stabilize the head in both the training and actual fMRI sessions. In multiple 256-s practice runs, subjects were trained to keep their head completely still while performing the simulated reach-to-eat task.

Image acquisition and preprocessing
Subjects were scanned with an 8-channel head coil in a General Electric 3 T scanner at UCSD Center for Functional MRI. Two functional scans were acquired in each subject using a single-shot echo-planar imaging (EPI) sequence with the following parameters: bandwidth ¼ 62.5 kHz; flip angle ¼ 90 ; TE ¼ 30 ms; TR ¼ 2000 ms; field of view ¼ 200 Â 200 mm; matrix ¼ 64 Â 64; voxel size ¼ 3.125 Â 3.125 Â 3.5 mm; 31 axial slices; 128 TR per volume after discarding 8 dummy TRs. Two field map scans were then acquired with the same orientation and dimensions as the functional scans for correcting the geometrical distortion. Finally, an alignment scan was acquired using the fast spoiled gradient-echo (FSPGR) sequence at the same orientation as the functional images with the following parameters: field of view ¼ 256 Â 256 mm; matrix ¼ 256 Â 256; voxel size ¼ 1 Â 1 Â 1.3 mm; 106 axial slices. In a different session, two additional sets of highresolution structural images (FSPGR; field of view ¼ 256 Â 256 mm; matrix ¼ 256 Â 256; voxel size ¼ 1 Â 1 Â 1 mm; 160-170 axial slices) were acquired for each subject.
The geometric distortion in the functional images was corrected using the field map scans and protocols provided by UCSD Center for Functional MRI (http://fmri.ucsd.edu/Howto/3T/fieldmap.html). Minor motion artifacts in the functional images were corrected using the 3dvolreg tool of the Analysis of Functional NeuroImages (AFNI) software package (Cox, 1996). The cortical surfaces of each subject were reconstructed from averaging two sets of high-resolution structural images using the FreeSurfer software package . Functional images were registered with the cortical surfaces by blink comparison between functional and structural images using an initial transformation matrix obtained by registering the alignment images to the high-resolution structural images using tkregister in csurf (https://mri. sdsu.edu/sereno/csurf).

Fourier-based analysis
To find task-related periodic brain activations, a 128-point discrete Fourier transform (DFT) is applied to the fMRI time series x m (t) of Voxel m by: where X m (ω) is the complex component at frequency ω in the power spectrum between 0 and 63 cycles per scan (cps), and jX m (ω)j and θ m (ω) are the amplitude and phase angle respectively. The stimulus frequency is defined as ω s (16 cps), and the remaining frequencies are defined as ω n (52 bins between 0 and 63 cps; excluding 0-2, 15-17, 31-33, and 47-49 cps). The "signal" and "noise" are defined as the complex components X m (ω) at ω s and ω n respectively. The F-ratio of this voxel can be obtained by computing the ratio between the normalized signal energy and noise energy (Dobie and Wilson, 1996): where df s ¼ 2 and df n ¼ 104 are the degrees of freedom of the signal and noise respectively. The p-value of F m is estimated by the cumulative distribution function F (2, 104) . In the present study, voxels containing strong periodic signals at the stimulus frequency (ω s ¼ 16 cps) with F (2, 104) > 7.39 (p < 0.001, uncorrected) were retained and their phase angles θ m (ω s ) were color-coded and rendered on the cortical surface using FreeSurfer (Fig. 2).

Independent component analysis
In the present study, a bite bar for immobilizing the head was not used because it would have blocked the hand-to-mouth movement. Consequently, the measured fMRI data contained periodic motion artifacts that were correlated with periodic reach-to-eat movements. To resolve this confound, we used independent component analysis (ICA) to decompose fMRI data into spatially independent components (ICs), which is a statistical method for solving the blind source separation problem (Beckmann, 2012;Bell and Sejnowski, 1995;Calhoun and Adalı, 2006;Hyv€ arinen, 1999;McKeown et al., 1998aMcKeown et al., , 1998bMcKeown et al., , 2003). An ICA-pruned fMRI dataset was reconstructed following the rejection of ICs containing periodic or aperiodic motion artifacts or random noise.
Let X 2 R PÂQ be the matrix of the observed time series, which is the linear mixing of K hidden (unknown) sources C 2 R KÂQ . The rows of C, c(1, :), …, c(K, :), are the independent sources (or independent components, ICs), which are independent non-Gaussian random variables. Hence, where A ¼ ða 1 ; …; a K Þ 2 R PÂK is the unknown mixing matrix with full column rank (so that P ! K).
The objective of ICA is to find an unmixing matrix W 2 R KÂP such that C ¼ WX is the best estimate of the hidden sources C by maximizing the non-Gaussianity or minimizing the mutual information between the sources using higher-order statistical moments. The unmixing matrix W is the pseudoinverse of A, i.e., W ¼Â þ . Hence, the sources C can be estimated by: The reconstructed signals can be represented as: To reject artifact-contaminated ICs l 1 , …, l r , their contributions to the observed signals can be removed by setting the l 1 th , …, l r th columns ofÂ to zero vectors, and an "ICA-pruned" dataset is then reconstructed by: ICA can be performed temporally or spatially on a dataset. In temporal ICA, each row of X is the observed time series, i.e., X is channel-by-time; whereas in spatial ICA, each column of X is the observed time series, i.e., X time-by-channel (Calhoun and Adalı, 2006). For analyzing fMRI data, each 4-D (3-D volume by time) dataset is first reshaped into a 2-D matrix X, where each column is the time series of a voxel and each row is obtained by reshaping each 3-D volume into a vector. Spatial ICA is then applied to decompose the matrix X intoX ¼ÂĈ (Eq. (5)), where each column ofÂ is the time series of an IC, and each row inĈ is the spatial component map (Beckmann, 2012;Calhoun and Adalı, 2006;McKeown et al., 1998b;Smolders et al., 2007). Following artifact rejection, an ICA-pruned fMRI dataset is reconstructed using Eq. (6).

Time-frequency analysis
Time-frequency analysis can be used to analyze the temporal characteristics of periodic signals and phase angles in fMRI data . For Voxel m, the fMRI time series x m ðtÞ is multiplied by a moving window wðt À τÞ centered at time τ. The resulting time series is then subject to the discrete Fourier transform: ½x m ðtÞwðt À τÞ expð À jωtÞ ¼ jX m ðω; τÞ jexp½jθ m ðω; τÞ (7) where the phase angle is defined as θ m ðω; τÞ. The signal-to-noise ratio (SNR) at the stimulus frequency ω s and at time τ can be obtained by: For Voxel m in an ICA-pruned fMRI datasetX obtained by Eq. (6), the time series of SNR,S m ðω s ; tÞ, and the time series of phase angles, θ m ðω s ; tÞ, can be obtained by Eqs. (7)-(9).

Circular statistics
The distribution of an angular dataset, θ 2 È θ q É Q q¼1 ¼ fθ 1 ; …; θ Q g, can be modeled using circular statistics. The angles are first converted into unit-length complex numbers using Euler's formula, z q ¼ expðjθ q Þ. The first sample moment of the complex number is obtained by: where 0 r 1 is the mean resultant length and Àπ θ < π is the mean phase angle of this dataset (Fisher, 1993;Mardia and Jupp, 2000).
To estimate the mean phase angle of a set of complex numbers with sample moment can be obtained by Grabska--Barwi nska et al., 2012;Levick and Thibos, 1982;Ringach et al., 2002): where θ and r are the mean phase angle and resultant length of the complex numbers.

Modeling the distribution of phase angles
If the phase angles in an angular dataset are drawn from the circular uniform distribution, they are distributed between -π and π in an equally likely manner and their mean is undefined. Therefore, the circular uniform distribution serves as the null model for testing against other alternative models (Fisher, 1993). The Rayleigh test can be used to determine whether a set of angles θ 2 È θ q É Q q¼1 ¼ fθ 1 ; …; θ Q g are distributed uniformly or drawn from a unimodal distribution (Fisher, 1993;Mardia and Jupp, 2000). Assuming that the angles in this angular dataset are uniformly distributed in the interval ½ À π; πÞ, the upper tail probabilities of Qr 2 can be approximated by (Mardia and Jupp, 2000): (2) were rendered in color on the cortical surface of a representative subject, as indicated by the color ring (0 -360 for a 16-s cycle). All phase angles were adjusted with respect to the mean phase angle of sROI OPA (see Table 1). Black dotted circles indicate regions containing periodic motion artifacts which should be disregarded. LH: left hemisphere.
where r is the mean resultant length of this dataset. A small p-value indicates that the angles are drawn from a unimodal distribution rather than uniformly distributed.
The 100ð1 À αÞ% confidence interval (CI) for the mean phase angle, θ, of an angular dataset θ 2 È θ q É Q q¼1 ¼ fθ 1 ; …; θ Q g can be estimated from bootstrapping the dataset (Fisher, 1993). Prior to bootstrapping, θ is The mean vector z * 0 and covariance matrix C 0 of fz q g Q q¼1 can be computed by: and Let V 0 be the square root of C 0 ; i.e., be the shifted angles of the corresponding vectors. To obtain z is formed by drawing the samples from fz q g Q q¼1 with replacement. The mean vector of fζ ðbÞ q g Q q¼1 , ζ ðbÞ , and the covariance matrix C b , are obtained from: and z * b can then be obtained by is rotated into the interval of ½ À π; πÞ when necessary.
The angles θ * 1 ; …; θ * B are then sorted in ascending order to yield a new set The 100ð1 À αÞ% confidence interval for θ can then be obtained by , where u is the integer part of 1 2 ðBα þ 1Þ (Fisher, 1993). The range of 100ð1 À αÞ%-CI for the mean phase angle can be computed by: (20) where B is the number of bootstrap samples. For Voxel m inX, the mean resultant length, r m , and mean phase angle, θ m , were obtained from its complex SNR time series,S m ðω s ; tÞ, using Eq. (11). The Rayleigh test was applied toS m ðω s ; tÞ to determine whether the phase angles are uniformly distributed using Eq. (12). If the phase angles from Voxel m are not uniformly distributed, the range of 95%-CI for the mean phase angle, ϑ ð95Þ m , of this voxel was assessed by bootstrapping the samples using Eqs. (13)-(20). Note that in the present study, z q in Eq.

Data processing pipeline
For each of 18 functional scans (9 subjects), initial maps of periodic activations were obtained by Fourier-based analysis (see Section 2.4 and the results of a representative scan in Fig. 2). To remove artifacts and random noise, spatial ICA was applied to each fMRI dataset using the FastICA algorithm (Hyv€ arinen, 1999;Hyv€ arinen and Oja, 2000), which was implemented in the FSL (FMRIB Software Library) MELODIC (Multivariate Exploratory Linear Optimized Decomposition into Independent Components) toolbox (Beckmann and Smith, 2004). Each 4-D fMRI dataset (64 Â 64 Â 31 Â 128) was reshaped into a 2-D matrix (128 Â 126976; where each column is the fMRI time series of a voxel), and then decomposed into a new dataset containing 30-54 spatial ICs (across 18 scans). For each ICA-decomposed dataset, ICs were manually labeled based on their activation maps and activation time series (Supplementary Figures S1-S5; McKeown et al., 1998b). ICs containing strong periodic or aperiodic motion artifacts or random noise were removed, and then an ICA-pruned datasetX was reconstructed with the original dimensions using Eq. (6) (Section 2.5).
Time-frequency analysis was applied to the time series of each voxel inX using a 32-point Chebyshev window (60 dB sidelobe attenuation; step size ¼ 1 point, total steps ¼ 97 points) in Eq. (7). A time series of the SNR,Sðω s ; tÞ, and a time series of the phase angle,θðω s ; tÞ, at the stimulus frequency were then obtained by Eqs. (8) and (9) (Section 2.6). The mean phase angle, θ, and resultant length, r, of the time seriesSðω s ; tÞ were computed using Eq. (11) (Section 2.7).
The Rayleigh test was used to determine whether the phase angles iñ θðω s ; tÞ of each voxel are uniformly distributed (Section 2.8; Eq. (12)). A voxel with p ! 0.001 (Bonferroni corrected) was rejected for containing uniformly distributed phase angles. For each retained voxel, the range of 95%-CI for the mean phase angle, ϑ ð95Þ , was then obtained by bootstrapping with 200 repetitions using Eqs. (13)-(20). If ϑ ð95Þ is larger than an empirically selected threshold (e.g., 12 for Subject 1, Scan 1; Fig. 3), θðω s ; tÞ is considered to have higher jitter in phase angles and this voxel will be rejected. Finally, the mean phase angle of each retained voxel was displayed in color on inflated cortical surfaces using FreeSurfer (Fig. 4).

Results
The objective and results of each analysis are illustrated in detail for Scan 1 of Subject 1 (Figs. 2-7). The results are shown only for the left hemisphere which is contralateral to the right visual hemifield and the right hand/arm involved in the simulated reach-to-eat task. Fig. 2 shows an initial map of periodic activations obtained by conventional Fourierbased analysis. This map demonstrates the periodic artifacts that cannot be separated by the conventional linear analysis. These artifacts were then identified and removed using ICA. Examples of task-related ICs and artifact ICs are shown in Figs. S1-S5. Following artifact rejection, time-frequency analysis and circular statistics were used to model the distribution of phase angles in periodic fMRI signals. Fig. 3 shows a map of the range of 95%-CI for the mean phase angle. Voxels were selected based on a mask obtained by setting an empirical threshold (12 ) in Fig. 3. Mean phase angles in the retained voxels from Fig. 3 are colorcoded and rendered on inflated and flattened cortical surfaces (Figs. 4  and S6). Fig. 5 shows keyframes of a movie (Video 1) created to unravel the complex spatiotemporal patterns of mean phase angles in Figs. 4 and S6. Fig. 6 shows 28 surface-based regions of interest (sROIs) outlined based on functional areas defined in our previous studies. The time courses of voxels enclosed in each sROI were averaged, and then sROIs were sorted by the mean phase angle of the average time course (Figs. 7 and S7). Table 1 summarizes the relative mean phase angle and delay for each sROI. Lastly, Figs. S8 and S9 show results of all single-subject scans, and Fig. 8 shows an inter-scan average map of mean phase angle.

Fourier-based analysis
An initial map obtained by Fourier-based analysis shows periodic activations forming major clusters in occipital, temporal, premotor, primary motor (MI), primary somatosensory (SI), secondary somatosensory (SII, including areas PV and S2; Disbrow et al., 2000;Huang and Sereno, 2007), and supplementary motor cortices (Fig. 2). However, periodic activations were also present in regions other than the visual and sensorimotor cortices. Some of the strong activations in orbital frontal lobe and anterior temporal lobe were motion artifacts resulting from periodic head movements (e.g., see the regions indicated by black dotted circles in Fig. 2). These artifacts were not corrected by the volume registration tool (AFNI 3dvolreg) but were later removed using ICA as detailed below.

Independent component analysis
Independent component analysis of fMRI datasets revealed ICs containing task-related periodic brain activations, periodic or aperiodic motion artifacts, or random noise. Task-related ICs were identified as containing clusters of periodic activations in primary sensorimotor, premotor, supplementary motor, lateral and medial posterior parietal, temporal, and/or occipital cortices (e.g., see task-related IC-10 and IC-22 of Subject 1, Scan 1 in Figs. S1 and S2).
Minor periodic head movements during a simulated reach-to-eat task resulted in periodic motion artifacts at the edges between the brain (bright) and non-brain background (dark) in 2-D functional images (data not shown). Surface-based maps of ICs containing these types of artifacts show periodic intensity changes in lateral prefrontal, orbitofrontal, anterior temporal, and lateral occipital regions (e.g., see artifact-related IC-1 and IC-20 of Subject 1, Scan 1 in Figs. S3 and S4; regions indicated by green circles). Other ICs containing aperiodic motion artifacts or random noise scattered within the brain were also identified based on their surface-based activation maps and activation time series. For each functional scan, all ICs containing the aforementioned artifacts or noise were rejected and an ICA-pruned fMRI dataset was reconstructed using Eq. (6).
In a direct-view setup where the subject's head is tilted forward, the ventral prefrontal and anterior temporal regions are located at the bottom of the scanned brain volume (axial slices). These regions are particularly sensitive to anterior-posterior (A/P) and pitch head movements. Furthermore, periodic activations in these regions are unlikely to result from sensorimotor related processes. Here, ICs containing periodic motion artifacts were further confirmed by comparing the phase angle of the 16-cps component of the IC activation time series with that of the A/P motion time series (output of AFNI 3dvolreg; see Figs. S1B, S2B, S3B, S4B, and S5). Head movements resulted in immediate intensity changes in the functional images with little to no delay in the activation time series, while a hemodynamic delay of several seconds was present in taskrelated brain activations. For example, the activation time series of an artifact component IC-1 and the A/P motion time series showed a small difference in their phase angles of the 16-cps component (Δθ ¼ 34.98 , equivalent to a delay of 1.55 s; Fig. S5). For another example, a smaller difference (Δθ* ¼ 21.2 , equivalent to a delay of 0.94 s) was present between an artifact component IC-20* (inverted waveform, see Fig. S5) and the A/P motion time series.

Circular statistics
For each voxel in an ICA-pruned fMRI dataset, a time series of phase angles was obtained by time-frequency analysis (Section 2.6) and then subject to the Rayleigh test (Section 2.8). The range of 95%-CI for the mean phase angle of each retained voxel is displayed in color on the cortical surface (Fig. 3). Voxels retained in the yellowish to reddish regions (with a range of 12 or smaller) are considered to contain stable periodic signals with low jitter in phase angles across time. These regions include the premotor cortex, MI/SI, SII, supplementary motor area (SMA), lateral occipital (LO) to middle temporal (MT) regions, areas V6/ V6A (Pitzalis et al., 2006(Pitzalis et al., , 2013, far peripheral representation of areas This single-subject map shows voxels with 0 -32 range of 95%-CI for the mean phase angle. Reddish to yellowish regions contain a range of 12 or smaller. V1 and V3, and ventral temporal-occipital cortex. The estimated mean phase angle θ of each retained voxel within the aforementioned regions is displayed in color on the cortical surface without spatial smoothing (Fig. 4). To find the relative temporal sequence among regions, all mean phase angles in the surface-based maps (Figs. 4-6) were adjusted with respect to the mean phase angle estimated from an sROI with the earliest visual response at trial onset (i.e., relative delay ¼ 0 s, relative phase angle ¼ 0 ; see Table 1).

Spatiotemporal dynamics across the cortical surface
To assist interpretation of a static map with color-coded mean phase angles (Fig. 4), a movie was created to show iso-phase contours traveling over the cortical surface for Scan 1 of Subject 1 (Video 1; Sereno and Huang, 2006;Windischberger et al., 2008). Keyframes of Video 1 were selected to illustrate sequential activations in occipital, temporal, parietal, and frontal cortices during a simulated reach-to-eat task (Fig. 5). To locate functional activations in the map and movie (Figs. 4 and 5; Video 1), group-average maps containing body-part representations (e.g., face, lips, and hands) and visual areas from the same subject pool in our previous studies were spherically morphed and then superimposed onto the cortical surfaces of Subject 1 (Fig. 6, left panel; Sereno, 2013, 2018;Huang et al., 2012Huang et al., , 2017. To illustrate brain regions with different activation delays, we outlined 28 sROIs for Scan 1 of Subject 1, each containing an approximately homogenous distribution of mean phase angles (Fig. 6, right panel and Table 1). The full names of abbreviations and detailed information of sROIs are summarized in Table 1. The time series of all voxels enclosed within each sROI were extracted and averaged (Fig. S7). The resulting average time series were then sorted by mean phase angles across sROIs and shown in a phase-sorted blood oxygen level-dependent (BOLD) image (Fig. 7A). A second BOLD image, created by reconstructing the 16-cps component from the power spectrum of the average time series of each sROI, clearly shows a traveling wave of periodic activations from the first to the last sROI (Fig. 7B).
The activation sequence of the 28 phase-sorted sROIs can be approximately divided into three stages based on temporal clustering (see color bar in Fig. S7): (1) sROIs 1-5 (early phase, reddish to yellowish); (2) sROIs 6-23 (intermediate phase, yellowish to greenish to cyan); and (3) sROIs 24-28 (late phase, bluish to purplish). The first stage involved visual areas (OPA, PPA, OTS, V7, and LO) that responded to the image swap at trial onset. The second stage began with activations in the superior parietal lobule (SPL1; Konen and Kastner, 2008) and dorsal premotor cortex (PMd) important for the planning of reaching, followed by areas involved in the execution of reaching (arm/hand representations in MI/SI, SMA, V6A, and VIPþ), and finally followed by areas involved in grasping (AIP, finger and lower face representations in SI, and ventral premotor cortex [PMv]). Furthermore, arm and hand movements near the lower right face (lower right visual field) activated visual motion sensitive areas CSv (Wall and Smith, 2008), MT, MST, VIPþ (parietal face area; Huang et al., 2017;Sereno and Huang, 2006), and far-peripheral visual field representations in areas V1, V3v, V3d, and V6. The third stage involved feeding related areas, including the human polysensory zone (PZ; Huang and Sereno, 2007) and face, lip, and tongue representations in MI/SI.
It is important to note that the discrete sROIs and stages discussed above were created to assist interpretation of the complex spatiotemporal brain dynamics continuously distributed across the cortical surface. Video 1 revealed several fine spatiotemporal patterns of traveling waves (as indicated by yellow arrows in Fig. 5), which could not have been observed from the temporal sequence of sROIs. First, at around 3.2 s, a wave emerged from area CSv and traveled anteriorly into SMA. Second, at around 3.6 s, a wave emerged from the arm representations in MI/SI, and then split into two waves dorsally and ventrally. The dorsal wave traveled into the upper arm and shoulder representations, which were involved in posture control during reaching movements. The ventral wave traveled into SI (finger representations) and AIP, which were involved in grasping movements. Third, at around 6.4 s, two waves emerged from the face/lip representations at the precentral and postcentral gyri, and then traveled inward to converge at the depth of the central sulcus (lip and tongue representations). Notably, the waves traveling progressively through hand, face/lip, and tongue representations in MI/SI revealed both the fine anatomical organization of hand and orofacial representations (see Fig. 6, left panel; Huang and Sereno, 2018) as well as the propagation directions of functional activations during feeding movements.

Inter-scan variation and average
Inter-scan variation in task compliance was assessed using the A/P motion time series, one of the motion parameters estimated by the AFNI Fig. 4. A map of mean phase angle estimated by circular statistics. This single-subject map shows the mean phase angles of retained voxels with a 0 -12 range of 95%-CI (Fig. 3). The color ring and offset of phase angles are the same as those in Fig. 2. See Video 1 for a traveling-wave movie created based on this map.  3dvolreg tool (Fig. S8). All 18 scans showed periodic head movements with different motion profiles (see average waveforms in the second and fourth columns in Fig. S8). Within each subject, the phase angles the of 16-cps component reconstructed from the power spectrum of the A/P motion time series showed little difference between the two scans (Δθ ¼ 19.6 AE 12.6 ; equivalent to 0.87 s AE 0.56 s). Within each scan, the A/P motion time series showed high periodicity as modeled by a cosine wave reconstructed from the 16-cps component. A nearly pure cosine wave of head movements suggested that the subject was performing selfpaced sequential movements at a consistent speed from trial to trial, e.g., see both scans of Subject 5 in Fig. S8. The spatiotemporal activation patterns in the maps of mean phase angle were somewhat variable across 18 single-subject scans (Fig. S9). To find activation patterns that were spatially and temporally aligned across scans, spherical morphing and surface-based complex averaging methods were used to average all single-subject maps of mean phase angle Hagler et al., 2007). For each scan, the surface-based map of mean phase angle (Fig. 4) was morphed into the spherical coordinates to register with an average sphere and then sampled to it using mri_surf2surf in FreeSurfer. For each vertex on the average sphere, an average complex number was obtained by vertex-wise vector averaging rexpðjθÞ (Eq. (11)) across 18 single-subject maps. The resulting spherical average map was then back-sampled onto the cortical surfaces of Subject 1 (Figs. 8 and S6; Video 2). Surviving areas in the inter-scan spherical average map include PMv, PZ, PMd, MI/SI (arm, hand, face, and lip representations), SII (PV/S2), 7b, AIP, VIPþ, PBA, SPL1, LIPþ, V7, MST, MT, LO, OPA, V3d/V6, V6A, CSv, SMA, and pre-SMA. The temporal sequence of activations among these areas was largely consistent with that in Subject 1, Scan 1 (Videos 1 and 2; Fig. S6). For example, the mean phase angles of SPL1, CSv, V6A, and VIPþ were ahead of those of PMd, SMA, pre-SMA, and arm/hand representations in MI/SI (see Discussion).

Discussion
Eating with the hand and arm is a trivial task for children and adults at all ages (Flindall and Gonzalez, 2013;Foroud and Whishaw, 2012;Quinlan and Culham, 2015;Sacrey et al., 2012;Sveistrup et al., 2008). Some populations with neurological disorders or brain injuries exhibit impaired movements while eating (Doan et al., 2008;Foroud and Whishaw, 2010;Hung et al., 2012). Due to technical constraints and  (Fig. 4) overlaid with contours and labels of visual and sensorimotor areas mapped in previous studies Sereno, 2013, 2018;Huang et al., 2012). PBA: parietal body area (including shoulder and leg representations). See text and Table 1 for the other abbreviations. Right panel: The same map overlaid with contours and indices of selected sROIs, each of which was manually outlined to enclose a cortical patch containing an approximately homogenous distribution of mean phase angles. The indices of phase-sorted sROIs indicate the sequence of activations, from the earliest (#1) to the latest (#28) in Figs. 7 and S7.
limited space in the MRI scanner, the reach-to-eat task has rarely been studied with fMRI in healthy or clinical populations. In the present study, we demonstrated experimental design and setup for investigating the spatiotemporal brain dynamics during a simulated reach-to-eat task in healthy adults. Head motion artifacts are inevitable for sensorimotor fMRI experiments in all populations. Across subjects in the present study, the time series of six motion parameters (output of AFNI 3dvolreg) showed minor periodic displacement in the superior-inferior (S/I) and anterior-posterior (A/P) directions, which were temporally correlated with periodic reach-to-eat movements. There was little to no displacement in the left-right (L/R) direction, suggesting the head was well stabilized laterally by foam padding. Minor head movements in the S/I and A/P directions were unavoidable because eating food in a realistic setting involves coordination among the head, trunk, arm, and hand (Foroud and Whishaw, 2012;Hung et al., 2012;Sveistrup et al., 2008;van der Kamp and Steenbergen, 1999). Although subjects had been instructed to keep   their head still, it was not possible to completely prevent the head from involuntary leaning forward to align the mouth with the incoming food.
While not preferable, allowing minor but correctable motion artifacts makes it possible to investigate the neural mechanism of head advancement while eating food. Furthermore, the task-related motion artifacts also provide additional information for functional data analysis (see discussion immediately below). Motion artifacts during the simulated reach-to-eat task were corrected in two stages. First, volume registration (AFNI 3dvolreg) yielded the time series of estimated motion in six degrees of freedom, which revealed not only the linear trends but also periodic fluctuations of head movements. Across subjects, these time series showed gradual linear drifts without abrupt motion, suggesting that all subjects were able to remain immobile for most of the experiment. However, volume registration was unable to completely correct and remove periodic motion artifacts correlated with periodic reach-to-eat movements. Second, independent component analysis was used to identify ICs containing task-related periodic brain activations, ICs containing periodic or aperiodic motion artifacts, and ICs containing scattered random noise. Periodic motion artifacts usually occurred at the edges between the brain and non-brain regions in functional images, which were likely to result from minor head displacement in the anterior-posterior direction and pitch head movements. These movements introduced immediate changes in image intensity at brain edges, as evident in the activation time series of ICs containing periodic motion artifacts (Fig. S5). In the present study, we simply removed ICs containing motion artifacts and noise. For ICs containing task-related periodic brain activations distributed across the cortical surface, further studies are required to analyze and interpret the spatiotemporal patterns of each IC in detail.
Previous event-related fMRI studies have used GLM-based designs and analyses to map brain regions that are involved in visually-guided actions, such as reaching, grasping, and pointing (e.g., Cavina-Pratesi et al., 2010Gallivan et al., 2011Gallivan et al., , 2016Rossit et al., 2013). The typical duration of each action event in these studies is about 2-4 s. Event-related activation in a brain region can be characterized by a hemodynamic response function (HRF), which reveals when the activation rises, peaks, falls, undershoots, and returns to the baseline. However, conventional methods using a time-fixed HRF model typically compare the peak magnitude rather than the latency of HRF across different conditions (actions) and brain regions. While time-resolved fMRI studies have used time-shifted HRF models to reveal the different activation latencies among brain regions (Formisano et al., 2002;Liao et al., 2002;Weilke et al., 2001), the time delays used in these models are usually limited to the multiples of an integral TR (e.g., 2 s, 4 s, 6 s, and so on). Although it is possible to create time-shifted HRF models with non-integral delays at high temporal resolutions, both the models and fMRI time series need to be up-sampled accordingly (e.g., from 2 s to 0.1 s per sample).
In the present study, we used the phase-encoded design and Fourierbased analysis to unravel sequential activations among brain regions during a simulated reach-to-eat task. The Fourier transform of fMRI time series revealed the phase angles (or delays) of periodic activations at 16 cps cycles per scan (16 s per cycle). Although the 16-s cycle was shorter than typically used in a phase-encoded mapping experiment (e.g., 32 s or 64 s), the resulting phase angles were continuously distributed around the full cycle (0 -360 ). This is because the intensity of functional images was recorded with a 16-bit quantization level, which allows timefrequency analysis and circular statistics to estimate phase angles (time delays) of periodic signals at a temporal resolution (e.g., 1 or smaller) much higher than the typical sampling rate of a phase-encoded fMRI time series (e.g., a 2-s TR in a 16-s cycle is equivalent to 45 in a 360 cycle). Furthermore, the 16-s cycle is sufficient for the subject to slowly perform a simulated, self-paced, and continuous reach-to-eat task, following a single visual cue (image swap) of action onset.
One of the disadvantages of a continuous and self-paced task is that the timing of movements could vary slightly from trial to trial. Although we could not record the subject's hand and head movements behind the direct-view screen, we were able to analyze the periodicity and regularity of task-related head movements using the time series of motion parameters generated by the AFNI 3dvolreg tool (Fig. S8). Across scans and subjects, the A/P motion time series showed strong periodic signals at the stimulus frequency (16 cps), suggesting that the timing of head movements was consistent across trials within each scan. Alternatively, it is possible to present an explicit visual or auditory cue for each action component (e.g., looking, reaching, grasping, or feeding). However, it is unnatural to perform a discontinuous reach-to-eat task with frequent pauses in a longer period (see future directions discussed below).
Conventional linear systems analysis (GLM or Fourier-based method) typically estimates the overall fit between the design matrix (square or sinusoidal waves) and the entire time series of each voxel (Boynton et al., Fig. 8. An inter-scan average map of mean phase angle. The results of spherical average of 18 single-subject maps of mean phase angle were rendered on the cortical surfaces of Subject 1. Colored regions indicate significant interscan activations (one-sample two-tailed t-test, t (17) > 3.97, p < 0.001, Bonferroni corrected for 18 scans). See Video 2 for a traveling-wave movie created based on this map. 2012). In our previous and present studies, we used time-frequency analysis to analyze the temporal fluctuation in periodic fMRI signals . Instead of a single F-statistic value that assessed the overall signal-to-noise ratio in the power spectrum (Sereno et al., 1995), time-frequency analysis yielded a time series of SNR and a time series of phase angles at the stimulus (task) frequency. Circular statistics was used to assess the stability of periodic fMRI signals based on the distribution of phase angles across time and estimate the mean phase angle of each voxel. A voxel was retained if it contained low jitter in phase angles across time, as assessed by the range of 95% confidence interval for the mean phase angle. The resulting map of mean phase angle showed periodic activations distributed across occipital, temporal, posterior parietal, primary and secondary somatosensory, primary motor, premotor, and supplementary motor cortices (Fig. 4). The overall extent of these activations was consistent with the regions activated by executed reaching and grasping, combined with the regions activated by passive tactile stimulation to the arm, hand, face, and lips in previous studies Filimon et al., 2007Filimon et al., , 2009Sereno, 2007, 2018;Huang et al., 2012Huang et al., , 2017Meier et al., 2008;Zlatkina et al., 2016). While previous studies have only shown static statistical maps of these sensorimotor representations, the present study further revealed the spatiotemporal dynamics among them.
The complex spatiotemporal brain dynamics during a simulated reach-to-eat task are best illustrated using a surface-based brain activity movie showing animated iso-phase contours (Video 1; Fig. 5). This approach is similar to tracking the progression of rainfall over time and space using an interactive Doppler radar map (e.g., https://weather. com/maps/usdopplerradar). Surface-based traveling waves facilitate understanding the hemodynamic response as a spatiotemporal process, rather than merely time-shifted HRF models of activations in spatially distinct regions of interest (Aquino et al., 2012;Lacy et al., 2016). For example, Video 1 showed a continuous wave of activations (or a gradient of mean phase angles) in a lateral occipitotemporal region including OPA, LO, and MT, and another wave in a cingulate region including CSv and SMA. Both the movie and average time series in phase-sorted sROIs (Figs. 6, 7, S7) revealed the temporal sequence of activations among brain regions involved in different stages of a simulated reach-to-eat task: (1) detecting the image swap activated scene/object selective areas (LO, OPA, PPA, and OTS) and V7 (reddish to yellowish); (2) planning and executing reaching movements toward the food image activated CSv, SPL1, PMd, c-SMA, V6A, VIPþ, MT, and arm representations in MI/SI (yellowish to greenish); (3) grasping the virtual food activated AIP, far-peripheral representations in V1, V3v, V3d and V6, r-SMA, hand (finger) representations in MI/SI, face representations in SI/SII, and MST (greenish to cyan); and (4) hand-to-mouth (feeding) movements activated PMv, PZ, and face, lip, and tongue representations in MI/SI (bluish to purplish). While it is possible to sort the phase angles of sROIs within each stage, sROIs in close temporal proximity to each other suggested that they could be functionally connected and coherent. For example, areas CSv, SPL1, and PMd were activated at approximately the same time (Video 1: frames 2.8-3.2 s), followed by simultaneous activations in MI/SI arm representations, MT, c-SMA, V6A, and VIPþ (see sROIs 6-13 in Figs. 6 and S7; Video 1: frames 3.6-4 s). The mean phase angles of reaching-related areas (sROIs 6-13) clearly preceded those of areas co-activated during grasping movements (sROIs 14-22; Video 1: frames 4.4-5.6 s).
While inter-scan variations were present in individual maps of mean phase angle, the spherical average map showed activation patterns that were spatially and temporally aligned across scans (Fig. 8). The coarse and fine spatiotemporal patterns of activations in the inter-scan average map were largely consistent with those in the representative scan discussed above (cf. Videos 1 and 2; Figs. 4, 8, S6). Among the surviving areas in the average map, SPL1, CSv, V6A, and VIPþ were activated the earliest, followed by PMd, SMA, pre-SMA, and MI/SI (arm representations), and finally followed by PMv, AIP, face representations in S2 and 7b, and MI/SI (hand, face, and lip representations). Previous time-resolved fMRI studies have consistently shown activation in SMA and pre-SMA during the planning phase of simple finger movements, followed by activation in MI during the execution phase (e.g., Cunnington et al., 2005;Formisano et al., 2002;Richter et al., 1997;Weilke et al., 2001;Wildgruber et al., 1997;Windischberger et al., 2008). However, our results showed nearly simultaneous activations in SMA, pre-SMA, and MI (arm representations). This is likely because SMA and pre-SMA were involved in online planning and execution of slow arm transport, but less in grasping and feeding movements. Compared with tasks that involve discrete stages of planning and execution, a continuous reach-to-eat task reveals a more comprehensive picture of real-time spatiotemporal dynamics across the entire cortical surface (see future directions discussed immediately below).
This study is our initial step to explore the challenges in conducting a simulated, continuous reach-to-eat task in the MRI scanner and in analyzing fMRI data across the time, frequency, and spatial domains. In future studies, these methods could be expanded in the following directions. First, real food such as dried nuts could be used to provide somatosensory and proprioceptive feedback to the fingers, mouth, and tongue during the task (Freud et al., 2018). However, subjects may only touch their lips with the food without actually eating it in order to avoid chewing and swallowing induced motion artifacts as well as choking risks. Second, MR-compatible video cameras could be used to record head, arm, and hand movements in the MRI scanner. Automatic video segmentation algorithms can then be used to analyze the exact timing of these movements. Third, future studies will also be able to take advantage of real-time video-based prospective motion correction during the functional scan. Fourth, explicit visual and/or auditory cues could be introduced in a discontinuous reach-to-eat task to study the cortical representations of different action stages (planning and execution) and components (reaching, grasping, and feeding). The exact timing of external cues and movement onsets/offsets recorded on videos would allow more accurate analysis of the temporal sequence of representations associated with each action stage and component. However, it is important to note that the complex spatiotemporal activation patterns observed during a self-paced and continuous reach-to-eat task may not be the linear summation (both spatially and temporally) of individual action representations mapped using a discontinuous task. Fifth, the phased-encoded fMRI data could also be analyzed using time-resolved multivariate (multivoxel) pattern analysis (Ariani et al., 2018;Gallivan and Culham, 2015;Gallivan et al., 2011Gallivan et al., , 2016. Lastly, time-frequency analysis and circular statistics could be applied to all frequencies in the power spectrum to estimate activation delays between brain regions (Sun et al., 2005).

Conclusions
It is challenging to study the reach-to-eat task with fMRI because of ergonomic constraints in the MRI scanner and task-related motion artifacts. In this study, we demonstrated the use of a direct-view setup for performing a simulated reach-to-eat task in a phase-encoded fMRI experiment. The subject periodically executed a compound action consisting of a continuous sequence of looking, reaching, grasping, and feeding within each 16-s cycle. These movements induced periodic motion artifacts in the functional images, which were identified and removed by independent component analysis. Time-frequency analysis and circular statistics were used to estimate mean phase angles of periodic fMRI signals. By sorting mean phase angles across surface-based regions of interest, we were able to show the real-time spatiotemporal brain dynamics among visual, multisensory, sensorimotor, and supplementary motor areas, which were involved in different stages of the reach-to-eat task. A movie was created to visualize sequential activations as continuous traveling waves over the cortical surface, which revealed both coarse and fine gradients of propagation as well as the interaction and coherence among brain regions. These complex dynamic patterns provided additional spatiotemporal information not obtainable by combining static cortical representations of action components individually mapped using GLM-based methods. In conclusion, this study demonstrates that phase-encoded design and analysis can be used to unravel sequential activations among brain regions at a temporal resolution higher than the native sampling rate of fMRI. The sequence of activations showed a surprisingly coherent traveling wave of activity similar to those seen in invasive experiments in the barrel cortex of awake rodents during whisking movements (e.g., Moldakarimov et al., 2018) as well as with voltage sensitive dyes in tangential slices (e.g., Wu et al., 2008). Although the dynamics of those activations are much faster than what we have observed, the similarities in spatial coherence are an intriguing avenue for future research.