Modelling Theta-Band Connectivity Between Occipital and Frontal Lobes: A Methodological MEG Study

Measuring functional connectivity between cortical regions of the human brain has become an important area of research. Modern theory suggests that brain networks exhibit non-stationarity, constantly forming and reforming depending on task demands. A robust means of determining effective connectivity in the short-lived neural responses that occur in event related paradigms would allow the investigation of event related cortico-cortical dynamics. We present such a mathematical model of wave propagation, motivated by current neuroscience literature, and demonstrate the utility of the method in a clinical sample of schizophrenia patients. MEG data were acquired in 10 patients with schizophrenia and 12 healthy controls during a relevance modulation task. Data were filtered into the theta band (4-8Hz) and source localised using a beamformer. The model was implemented using Fourier analysis methods which uncovered an event related travelling wave moving from the visual to frontal cortices. The model was validated using Monte Carlo phase randomisation and compared with another plausible model of wave propagation in the cortex. Results from the clinical sample revealed that wave speed was modulated by task condition and patients were found to have less difference between conditions (ANOVA revealing a significant interaction between group and condition, p<0.05). In conclusion, our method provides a simple and robust means to investigate event related cortico-cortical brain dynamics in individuals and groups in the task positive state.

number of robust large-scale distributed networks of connectivity; some associated with sensory systems (e.g. the visual and sensori-motor networks) and others associated with attention and cognition (e.g. the default mode and dorsal attention networks).
Using electrophysiological measurements, several authors have proposed that phase synchronisation of oscillatory activity serves as an effective means of cortico-cortical communication Singer and Gray, 1995) in both short range (Eckhorn et al., 1988;Engel et al., 1991; and long range connections Fries, 2005;Hipp et al., 2012;Salinas and Sejnowski, 2001;Varela et al., 2001). Though useful, current methods primarily aim to identify aggregate brain networks observable over large time windows, in many cases lasting several minutes. However, human brain function relies on our ability to react to stimuli on a rapid timescale. This must necessitate the rapid formation and disintegration of transient brain networks, defined by short lived connections that are superimposed upon the distributed networks observable using traditional static metrics (Engel et al., 2013). fMRI is an indirect measure of brain activity based upon slow haemodynamics and thus is not well suited to the investigation of these transient brain states.
We aimed to use the relatively short-lived induced envelope from a seed location, and predict activity at other brain regions using a transfer function. This transfer function takes into consideration the non-linear effects of temporal delay, thereby capturing the transient effective connectivity associated with a particular cognitive task. The model presented here is informed by relevant knowledge of neuronal functioning showing that oscillatory coupling between regions results in effective connectivity and cortical communication.
Communication between cortical regions may be achieved by synchronised firing of neuronal populations, since the integration of such coincident firing at the target population increases the likelihood of a neuronal response (Abeles, 1982;König et al., 1996;Salinas and Sejnowski, 2001). Therefore, tightly synchronised oscillations are more effective than non-synchronised outputs at transferring information from one region of the cortex to another (Alonso et al., 1996;Azouz and Gray, 2000;Bruno and Sakmann, 2006). Engel et al. (2013) have described short lived spatiotemporal patterns of connectivity as intrinsic coupling modes (ICMs) categorised loosely by two mechanisms of interaction; the first mechanism involves a linear fixed phase relationship between the oscillations generated by spatially distinct brain areas (Gow Jr et al., 2008;Gross et al., 2001;Ioannides et al., 2000;Jerbi et al., 2007;Nolte et al., 2004;Schlögl and Supp, 2006;Schoffelen and Gross, 2009;Tass et al., 1998). The second involves temporal correlation between the amplitude envelopes of neural oscillations in different regions (Brookes et al., 2011a(Brookes et al., , 2012Bruns et al., 2000;Hall et al., 2013;Hipp et al., 2012;Liu et al., 2010;de Pasquale et al., 2010). Both forms of coupling identify networks of functional connectivity in the human brain. Furthermore, both forms have a degree of spatial concordance with the networks observed using fMRI (Brookes et al., 2011a;Marzetti et al., 2013;de Pasquale et al., 2010;Schölvinck et al., 2010).
However, the precise mechanism by which amplitude correlation occurs is not known. One mechanistic interpretation is that state changes in distributed cortical networks are modulated on slow timescales by neuro-modulatory input (Leopold et al., 2003;Munk et al., 1996). However, on shorter timescales, i.e. during a cognitive task containing short trials, amplitude modulation is likely to be driven by transient synchronisation of local cortical oscillations (Neuper et al., 2006;Pfurtscheller and Aranibar, 1977;Pfurtscheller and Lopes da Silva, 1999). The true picture is likely a combination of the two. However, for building our model, we will assume that transient connections between cortical regions during a task are driven by inter-cortical communication. When two or more neuronal populations are engaged in coherent activity, phase differences between sending and receiving populations may regulate the flow of information (the effective connectivity) (Fries, 2005;Womelsdorf et al., 2007) a proposition that forms the basis of effective connectivity measures in the EEG and MEG such as dynamic causal modelling (DCM) (Friston et al., 2003;Stephan et al., 2007).
In the stimulus induced response this directed form of connectivity may also manifest in the envelope as correlation of event related synchronisation (ERS), even in the absence of long range phase coherence (Bruns et al., 2000). Therefore, inter-cortical delays in transient ERS and event related desynchronisation (ERD) might provide a phase independent marker for effective connectivity.
Importantly, this theory implies that greater synchronisation corresponds to communication that is more efficient. Locally, this should be evident in the envelope as more rapid ERD and ERS in the post stimulus response (in other words a more rapidly modulated envelope in both communicating regions). In terms of cortico-cortical communication, this increased efficiency should result in a decrease in delay. Therefore, we make the prediction that temporal delay between cortical sites is proportional to the modulation frequency of the envelope.
In order to test a model based on these assumptions it is necessary to adopt a paradigm known to elicit long-range communication, and to elicit reliable components in the stimulus-induced response. The theta band has been implicated in communication between the occipital and medial frontal regions during working memory (Sarnthein et al., 1998) and visual attention paradigms (Gregoriou et al., 2009), suggesting a role for theta in long range communication between distant cortical sites.
Furthermore, the amplitude of post stimulus theta ERS scales with task difficulty (Brookes et al., 2011a(Brookes et al., , 2011bJensen and Tesche, 2002) supporting the idea that theta ERS is somehow functionally significant in visual working memory. Anatomically, the visual processing pathway is comprised of a cortical hierarchy that sub-serves cognitive function by participating in both feed-forward and feedback cortical connections. Feed-forward signals convey sensory information, while feedback mechanisms provide modulatory input (Bastos et al., 2015a(Bastos et al., , 2015bFelleman and Van Essen, 1991;Markov et al., 2014). Furthermore, sensory-guided goal-directed behaviour necessary for the completion of most cognitive tasks, relies on both feed-forward (Brookes et al., 2011a;Donner et al., 2009;Glimcher, 2003;Mazaheri et al., 2009Mazaheri et al., , 2010 and feedback networks (Donner et al., 2008(Donner et al., , 2009Driver et al., 2009;Gold and Shadlen, 2001;Horwitz et al., 2004;Kastner and Ungerleider, 2000;Kim and Shadlen, 1999;Moore and Armstrong, 2003;Nienborg and Cumming, 2009;Ress and Heeger, 2003;Serences and Yantis, 2006;Zanto et al., 2011). Bastos et al. (Bastos et al., 2015a) have recently shown that, in the primate visual system, feed-forward and feedback influences are carried by distinct frequency ranges of neural oscillations. Specifically, theta (~4Hz) and high gamma band (~60-80Hz) oscillations were involved in routing of information from V1 to higher processing areas, while feedback influences were carried by beta (~14-18Hz) oscillations. Inter-regional delay of power fluctuations between several functional networks have also been observed during a visual working memory study in humans (Brookes et al., 2011a;Liddle et al., 2016) suggesting that the relevance modulation task is a useful paradigm for investigating striatal-frontal connectivity. It was found, in both of these studies, that trial averaged induced responses revealed a burst of theta oscillations peaking at ~200ms post stimulus onset in the visual cortex. A similar theta waveform was observed in the frontal regions, but was delayed relative to the stimulus by ~400ms. If this observable theta burst is indeed representing a feed-forward influence within a hierarchical network it follows that a direct correlation should be observable between the theta waveform within these separate regions, if the mapping function between those two regions can be determined accurately.
Here, we aim to investigate the precise relationship between visual and frontal theta-band envelopes.
Using MEG, we measure the theta band waveforms that are induced across multiple brain locations by a cognitive stimulus. We formulate two simple models of propagation of the theta envelope waveform from the occipital cortices to anterior brain regions. In the first model, we assume that wave components at all frequencies travel through the brain with the same velocity, such that the theta waveforms measured in the frontal regions are simply time-delayed reconstructions of the equivalent theta waveform in the visual cortex. In the second model, we employ a model to test the prediction that higher frequency envelope modulations will travel through the brain with a shorter delay. As a result, the envelope modulation frequency will be linearly related to the propagation delay, leading to a temporally distorted theta waveform in distant cortical regions. Note that this is not the same as measuring phase lag in the carrier signal. We present evidence that the latter model offers a useful means to investigate connectivity between distal brain regions (specifically visual to frontal). We also show that the parameters characterising the model are affected by the cognitive task undertaken by the subject. Finally, we show that these same parameters differ in patients with schizophrenia compared to healthy controls.

Theory:
In order to assess connectivity between two distal brain areas we first begin with a simple description of the two models to be tested. Consider first a simple sinusoidal waveform, induced by an external stimulus (e.g. visual stimulation) within a neural network at brain location . This waveform is subsequently "transmitted" via some means to a spatially separate brain location, (see Figure 1A for an example).

If
, represents the measured waveform at location then taking the simplest possible waveform, we can assume , sin . [1] Note that amplitude, , and frequency, , are determined by the external stimulus driving the waveform (i.e. if the stimulus were a flashing light, then would likely be defined by the frequency of the flash) and the localised cortical response to that stimulus at brain location . Now assume that , is the transmitted waveform measured at brain location . Note that the signal has, via some unknown path, been transmitted between and over a distance , with some effective velocity . Mathematically we guess that Here, is the time taken for the wave to travel from to and so . [3] Therefore, , sin − . [4] Note here that unlike the general case for a physical wave (e.g. light waves), the parameters ω and are unrelated. ω is likely to depend on the stimulus and response function of brain area whereas represents conduction velocity in the brain.
The above treatment only accounts for a pure wave of a single frequency component. To generalise this, consider what happens if multiple different frequency components originate at . Here, , sin . [5] Again, assuming that a linear shift in the waveform between and we could write that , sin − . [6] Note that , the amplitude of each frequency component is allowed to vary with frequency. Recall also that for our simple case, and are constant. We now rewrite the expression for as In words we state that is a phase term which is directly proportional to frequency ω. Implicit in the model is that we assume that conduction velocity, , is frequency independent, meaning that the measured waveform at is simply a time shifted version of that measured at . This simple model is depicted in Figure 1B and, since is a first order polynomial in ω, we term this a first order phase shift model.
The assumption inherent in the first order phase shift model is that all amplitude modulation frequencies travel with the same velocity, . However, alternative models may be proposed where other factors affect propagation velocity. One possible factor is the synchrony of the cortical population. Here we assume that a rapid modulation of band-limited power reflects a rapid change in synchrony of the local neural population. Neurons within the local population entrain distant connected neurons resulting in a rapid change in power at the connected region. Therefore both regions will show rapid modulation of power with a time delay that is proportional to the modulation frequency. This effect is independent of axonal conduction velocity which is assumed here to be relatively fast and of constant velocity. On different trials, and within different neuronal subpopulations, amplitude modulations of varying frequency add linearly resulting in a more complex average waveform.
For a simple idealised envelope with sinusoidal time course: Where c is an unknown constant. Note that we do not suggest that physical conduction velocity through white matter changes with wave frequency, rather we use as an effective velocity which accounts for relative delay in action potential generation caused by different group delay. To denote this, we employ the term !! . Substituting Equation [9] into the above model we can write: Or simply: Where ' = /" and is a constant for all frequency components. In other words, the phase shift is independent of frequency (i.e. ' is a zero-order polynomial). We term this model a zero-order phase shift. In fact, the zero order phase shift results not in a time shifting of the original waveform as depicted in Figure 1B, but rather a distorted version of the original waveform. This depicted in Figure   1D which shows two flow charts detailing the first order and zero-order phase shift models of transmission between two brain regions.

Materials and Methods:
Participants and Paradigm: MEG data were acquired in 22 subjects. This comprised 10 patients satisfying the DSM-IV criteria for schizophrenia (Leckman et al., 1982) (Liddle et al., 2002), and 12 healthy controls not differing significantly from the patient group in mean age, sex or parental socioeconomic status (Pevalin and Rose, 2002). Nine patients (75%) were receiving anti-psychotic medication. Patients below the age of 18 or over 50, or with an IQ score of below 70 on the Quick Test (Ammons and Ammons, 1962), were not considered for the study. Healthy controls were recruited from the community via advertisements and were free from Axis I psychiatric disorders. The study was given ethical approval by the National Research Ethics Committee, Derbyshire, United Kingdom. All participants gave written informed consent.
All subjects took part in a "Relevance Modulation" task designed specifically to modulate the relevance to the subject of incoming visual stimuli. Note that this task has been employed in previously published studies (Brookes et al., 2011a) and we include its description here again for completeness.
The subject was shown a series of visual stimuli (duration=800ms) in which pictures of butterflies were alternated with pictures of ladybirds. Inter stimulus intervals (ISI's) were jittered, a single stimulus being presented, on average, every 2s. A single block comprised presentation of 40 alternating stimuli, followed by a resting phase. Blocks were 120s in total. At the beginning of each block, subjects were informed whether the relevant stimuli were the butterflies or the ladybirds. In blocks where butterflies were relevant, the subject's task was to respond when the butterfly on the screen matched an example butterfly (shown at the start of the block) on three dimensions (shape, outer wing colour; inner wing colour); in blocks in which the ladybirds were relevant (LB blocks), the subject's task was to respond if the number of red ladybirds in the picture was equal to the number of yellow ladybirds. The number of ladybirds ranged from 4 to 6, and they were presented in pseudo-random positions and orientations. In both block types, the "target" stimuli were rare (probability = 0.05), and these "target" trials were excluded from analyses so as to avoid confound from a motor response. Eight blocks were presented to each subject, in the order: BF, LB, LB, BF, LB, BF, BF, LB. Both task types (responding to a butterfly that matched the example; responding to equal numbers of red and yellow ladybirds) were designed to require close attention to the relevant stimulus.

Data Acquisition and Preprocessing
MEG data were acquired using the third order synthetic gradiometer configuration of a 275 channel MEG system (MISL, Coquitlam, Canada), at a sampling rate of 600Hz and using a 150Hz low pass anti-aliasing filter. Three head position indicator (HPI) coils were attached to the subject's head at the nasion, left preauricular and right preauricular points. These coils were energised with a subject inside the scanner to allow localisation of the head relative to the geometry of the MEG sensor array. The surface shape of the subject's head was then digitised using a 3D digitiser (Polhemus, Isotrack) relative to the coils. The surface shape of each subject's head and the coil locations were then coregistered with an anatomical MR image, acquired using a Philips 3T Achieva MR system at 1x1x1mm 3 isotropic resolution using an MPRAGE sequence.
In order to project sensor space MEG data into brain space a scalar beamformer spatial filtering approach was employed (Robinson and Vrba, 1999;Van Veen et al., 1997). Beamformer weighting parameters were computed on a regular 8mm grid spanning the whole brain space. A different covariance matrix was computed for each frequency band. The forward model was based upon a multiple local sphere head mode (Huang et al., 1999) and the dipole model described by (Sarvas, 1987). Multiplication of the theta band filtered MEG data by the beamformer weights generated a timecourse of electrical activity for each voxel. Theta-band envelope time-series were then generated Hilbert transformation. These envelopes were averaged across all of the relevant and irrelevant trials in the dataset in a time window spanning 0+ < < 2+ where represents peri-stimulus time.

FOPS and ZOPS Model Generation
Our primary aim, as detailed in our introduction, was to test whether the zero order phase shift (ZOPS) model was better able to measure transmission of a waveform through the brain than the first order phase shift (FOPS) model. To do this, we first take the measured waveform, , at brain region and act upon it using simple Fourier analysis to simulate the waveform at . The simulated waveform can then be compared to the measured waveform in order to test which of the two models best fits the data. We denote the Fourier transform of , as . , which can be written: Where, represents the amplitude and the phase of each frequency component in Fourier space. We now simulate the signal at based upon the signal from , using the FOPS model as: 2 34 , = . 5 / 0 1 51 6 [13] Note here . 5 represents the inverse Fourier transform, and we use the 'hat' notation 2 to represent the simulated version of . Note also that as required for the FOPS model, = 7 where p is a constant. Similarly it is possible to simulate the signal at based upon the signal from , using the ZOPS model as: Where ' is a constant, meaning that to generate the simulated ZOPS signal we modify the phases of each frequency component by the same amount, irrespective of frequency.

Model testing and statistical analyses
In order to compute the similarity between the model and the data, correlation coefficients were computed between the measured theta timecourse and the modelled theta timecourse 2 . In order to generate the model, a time series, , was taken from a seed location in primary visual cortex and used, alongside either Equation 13 or 14, to generate the modelled timecourses 2 .
For the zero order phase shift, ' was allowed to vary, taking 100 equally spaced values in the range :0 < < 2;<. A set of 100 correlation coefficients was calculated for every voxel in the brain and in all cases the highest correlation across the 100 phase values was extracted and used in further analysis (i.e. we select the phase shift value that best fits the data for any one single voxel).
For the first order phase shift, the phase change effectively amounts to a shift in time of the original seed time series. 100 iterations of the model time series were therefore computed with the time shift equally spaced between 0s and 2s. Again a set of 100 correlation coefficients was calculated, for every voxel in the brain, and the highest correlation across the 100 time shift values was extracted (i.e. in this case we select the time shift value that best fits the data for any one single voxel).
A potential limitation of the correlation approach is that for MEG, signals derived from multiple spatially separate voxels are not necessarily independent. This "signal leakage" between voxels is a direct result of the ill posed MEG inverse problem and is particularly problematic for voxels close to the centre of the head. In such cases the signal to noise ratio is very low and the projected timecourse can often represent a heterogeneous mix of the cortical signals. For this reason, in addition to correlation coefficients we also computed covariance. Computing the covariance between the model and the data not only accounts for similarity between and 2 but also the variance in . In beamforming, the variance of the signal is downweighted according to the norm of the weights vector (see Hall et al 2013) meaning that in deep regions, where leakage is high, the variance of 2 will be relatively low. We would therefore hypothesise that the covariance metric would generate a more accurate spatial map than the correlation metric.
To generate the correlation and covariance images, each subjects individual MRI was segmented and the brain extracted using the brain extraction tool (BET2) in FSL (Smith, 2002). The brain extracted MRI for each subject was then coregistered to MNI space using FLIRT (Jenkinson and Smith, 2001) and the same transform was applied to the functional (beamformer derived) data which were averaged across subjects in standard space. The covariance and correlation maps were then computed on an 8mm grid in MNI space using the subject averaged data. All covariance and correlation maps were interpolated to a 1mm 3 isotropic MNI grid, overlaid on to a standard MNI brain (MNI152 1mm brain, available in FSL) yielding 100 covariance and 100 correlation maps (i.e. one map was available for each of the 100 phase shifts (in the case of ZOPS) or time shifts (in the case of FOPS). Animations were generated by sequencing these images from low to high phase shift. Finally, maps of optimum phase for every voxel in the brain were generated.
Calculation of statistical confidence intervals for each method, and comparisons between methods, were approximated by Monte Carlo phase randomisation of real signals extracted from the control subjects. Firstly, timecourses were extracted using beamformer source localisation at 50 locations spanning the visual and frontal cortices, chosen to be representative of the whole head. For each iteration of the Monte Carlo simulation, the extracted signals were phase randomised in the Fourier domain using the same method as described in (Prichard and Theiler, 1994) to maintain the covariance structure between locations, and preserving zero phase lag correlations such as those caused by weights correlation (see Fig 1.02). Each surrogate dataset was analysed using the same methodology as that used in real data. The phase-randomisation process was repeated 1000 times for each location producing a null distribution of r values. All values from >90mm were grouped to produce a general null distribution. Figure 1.02. Results generated by phase randomisation simulations. The mean correlation of the null distribution is shown in blue, and the lower and upper 95% confidence limits are shown in grey below and above the mean. The effect of seed blur on the envelope correlation can be observed near 0mm. However, this flattens off at >30mm after which point the mean remains relatively constant. Note that FO method generates consistently higher mean correlation values under the null distribution than the ZO method. This difference needs to be taken into account when comparing procedures in real data.

Between Groups Comparisons (Clinical Data)
For the between groups comparisons of the controls and clinical sample, analyses were performed on individual subjects. Here the delay estimates for each model were calculated per voxel, and then averaged across distance bins producing 9 delay values (1 per distance bin) for each subject. These values were then checked for statistical normality. Our hypothesis was focused on the frontal cortex, so only the last two distance bins (encompassing the frontal cortices >120mm from the seed) were entered into an ANOVA model (2 x condition, 2 x group, 2 x distance).

Results:
Healthy Volunteer ZOPS Connectivity Analysis: Figure 2 shows results of our ZOPS connectivity analysis applied to data from healthy subjects, with a seed location in primary visual cortex. Analysis was limited to the averaged theta (4-8Hz) envelope timecourse. Figures 2A and B show axial and sagittal correlation images for phase values varying from 0 ο to 91 ο ; Figure 2A shows the case for the relevant stimuli whilst Figure 2B shows the case for irrelevant stimuli. (The locations of the axial and sagittal slices used are shown in Figure 2H.) It is noteworthy that the ZOPS model generates high correlation, with the voxels exhibiting the highest correlation appearing to move anterior, as phase value is increased. The peak correlation of r=0.97 was located in medial frontal cortex and p values (obtained using the Monte-Carlo phase randomisation method) survive a Bonferroni correction for all 4626 voxels tested (pcorrected<0.05). This was considered a very conservative significance test. The anterior movement of regions of high correlation between the data and the ZOPS model, as phase is incremented, generates the effect of a 'travelling wave' propagating through the head. Note that the correlation values are affected by task condition, with high correlation for relevant stimuli and lower correlation for irrelevant stimuli. Figure 2C and D shows the equivalent covariance images, derived by taking the correlation values in A-B and scaling by their normalised covariance. Note here that the images become more focal, with visual areas highlighted for low phase change and medial frontal cortices highlighted when the phase change approaches 90 ο . The increased focality of the highlighted regions likely results from the fact that variance normalisation, inherent in the correlation images, amplifies leakage close to the centre of the head; the fact that the variance of such regions is down-weighted by the beamformer means that they fail to appear when covariance is employed. Figures 2E and 2F contain example timecourse fits, for voxels located in the frontal cortex. The seed timecourse, Y , from the primary visual cortex, is shown in green, the test voxel timecourse, Y is shown in blue and the model fit, Y > , is shown in red. Note the excellent agreement between the model and the data, particularly for the relevant condition. Figure 2G shows phase maps, highlighting the best fitting phase angle for each voxel in the head. The relevant condition is shown on the left and irrelevant on the right. Note the higher phase values in the frontal regions for the relevant condition.
Note also the difference between the relevant and irrelevant conditions. Overall, the results presented in Figure 2 show good evidence that the ZOPS model generates good fits to measured theta-band envelope waveforms. Furthermore, the fitted phase angles change with location in the brain, with increased phase occurring in the frontal regions. Note also that the fitted parameters change with task condition.  Figure 4 shows a comparison of the zero-order and first-order phase shift methods, for both the relevant and irrelevant condition. Correlation values (derived from the optimum phase angle) are plotted against Euclidean distance between the seed and the test voxel. The ZOPS model is shown in blue, and the FOPS model in red. Figure 4A shows a single point for every voxel whilst in figure 4B, voxels are grouped into 9 equally spaced bins (representing Euclidean distance) and the mean and standard deviations across voxels in each bin displayed. The primary result shows that the ZOPS method is better able to explain variance in the data compared to the FOPS method for large Euclidean distances in the brain, specifically for those regions greater than ~120mm from the seed. As shown by the inset image, these voxels were primarily located in the frontal lobes. The difference between the two models, for the relevant condition and for voxels >120mm from the seed, was analysed statistically, revealing a significant difference in between methods (p<0.05). There was no significant difference between models for the irrelevant condition. The dotted line represents the 95 th percentile of the null distribution computed via the monte-carlo phase randomisation procedure (outlined in methods). Mean and standard deviation across voxels for each distance bin. The 95% confidence interval generated for each method is plotted as a dotted line. In both A and B, the relevant condition is shown on the left hand side and irrelevant on the right hand side.

Analysis in Schizophrenia and Correlations with Behaviour
Finally, we used the ZOPS model to look for interrelationships between the visual cortex and the frontal lobes in patients with schizophrenia. Figures 3A-D shows correlation and covariance maps, for the relevant and irrelevant conditions, in a group of 10 patients. Visual comparison between the control group in Figure 1 and the patient group in Figure 3 reveals differences between patients and controls in the magnitude of the values in the frontal cortex. Figures 3E and F show a comparison of the mean measured phase angle, plotted as a function of Euclidean distance from the seed, for the relevant (green) and irrelevant (red) conditions. Note that for group statistical comparisons, ZOPS analysis was performed on individual subjects giving a phase estimate per voxel, per subject. The results, illustrated in figure 3, show that in the control group, task condition had a large modulatory effect on delay. However, in the patient group delay did not appear to be modulated by task condition. For statistical comparison, the last 2 bins (those that encompass the frontal cortices) were entered into a repeated measures ANOVA [2 distance x 2 condition] with group as a between subjects factor, and the mean phase in each bin as the dependent variable. The analysis revealed a significant effect of condition (F=7.35, p<0.05), a significant effect of distance (F=5.09, p<0.05), and a significant interaction between group and condition (F=4.87, p<0.05). This demonstrates that overall, phase significantly changes with distance from seed location, and the presence of a task induces a significant change in phase in the controls but not the patients.
Finally, Figure 4 shows correlation between the measured phase from the ZOPS model, and behaviour; specifically we employed D-Prime as a measure of signal detection on the RM task (see methods). The results show that the degree of phase shift in the frontal cortex, during the relevant condition of the task, negatively correlates with discriminability (r=-0.48, p<0.05, two tailed).  included. The Y-axis shows D-Prime, a measure of the subject's ability to detect a target. Phase delay was negatively correlated with D-Prime (r=-0.48, p<0.05) supporting the hypothesis that task performance is related to inter-regional phase delay in the theta BLE.

Supplementary Figures
Supp. Fig. 1. The alpha band appears to show a 180 degree phase shift. i.e. an anti-correlation.
Supp Fig. 2. Beta band appears to show, a 90 degree phase lag between frontal and visual cortex. Here the beta-band power fluctuations in frontal cortex appear to precede those in visual cortex.

Discussion:
We have introduced a novel method to assess transient task induced connectivity between brain regions, based upon measurement of a phase delay between theta-band envelope signals measured across the entire cortex. We have compared two separate phase based models; in our first model, termed first order phase shift (FOPS), we assume that all frequency components in the theta envelope will travel through the brain at the same velocity, leading to a time shifted version of the seed envelope at the second brain region. In our second model, termed zero order phase shift (ZOPS), we propose that faster and more efficient processes will manifest in the envelope signal as higher frequency components. These components of the signal will travel with a higher effective velocity than the slower components, leading to a temporally distorted version of the seed envelope at distant brain regions.
We have shown that, at least for the cognitive task employed here, the ZOPS significantly better than the FOPS model, and that the ZOPS model performs significantly better than chance when using randomised surrogate data. Furthermore, we have shown that the fitted phase values are altered significantly by task condition, they are correlated with a behavioural measure, and the condition dependent modulation of phase seen in healthy subjects is not present in patients with schizophrenia.
This methodology represents a fundamentally new way to measure directed functional connectivity in the brain.

Correlation vs. Covariance
Correlation metrics of connectivity used here appeared to reveal that a large proportion of the brain is connected during the relevant condition of the task, with a spatio-temporal pattern similar to a travelling wave of theta-band activity. However, using covariance data in the theta-band allowed the discrimination of a task relevant network involving primarily the occipital lobes and medial frontal cortex. These findings appear to be consistent with previous research implicating theta oscillations in long range connectivity (Buzsaki and Draguhn, 2004;von Stein and Sarnthein, 2000;Varela et al., 2001) and research relating theta amplitude modulation to task performance in the frontal midline (Brookes et al., 2011b;Gevins et al., 1997;Jensen and Tesche, 2002;Luft et al., 2013). In the present study, mid frontal theta envelope variance is almost entirely explained (R 2 = 0.94 at the peak) by a ZOPS of the envelope from V1 indicating a close coupling of these two regions. We have also shown that inter-cortical latency is dependent on task condition, which had a significant modulatory effect on the dynamics of the theta envelope. The context dependent modulation of the theta band suggests some task dependent top-down modulation of feed-forward influences, in line with current predictive coding models of brain function and cognition (Arnal and Giraud, 2012;Friston, 2005). Abnormalities in the establishment of long-range connectivity via low frequency oscillations appear to be disrupted in patients

Neural Mechanisms Underpinning the ZOPS Model
We have used a mapping function with an occipital seed as the input to accurately model theta-band envelope time-courses across many distant cortical areas. However, it is unclear how the carrier signal of these envelopes supports this close relationship. It is known that connectivity in the envelope can exist without measurable phase coherence in the carrier signal, particularly in long range connections (Bruns et al., 2000). Furthermore, the effective connectivity observed here in the theta band was observed in the group average per condition. Phase coherence must necessarily be measured in the raw signal, to avoid cancellation caused by phase jitter. In future, we plan to use phase-based measures, such as granger causality to investigate the relationship between the carrier signal and the envelope.

Travelling Waves as Possible Mechanisms of Oscillatory Influence
Although the model put forward here has performed well in describing inter-regional variation in the band-limited envelope, the precise neurobiological mechanisms underpinning the apparent travelling wave observed are not fully understood. Inter-cortical phase coupling is widely thought to be the primary mechanism for communication between cortical sites Fries, 2005;Salinas and Sejnowski, 2001) and is modulated by stimulus context, task, or cognitive setting Siegel et al., 2012;Singer, 1999). However, our model is based on coupling of signal envelopes rather than oscillatory phase, and as discussed, the two are not directly related. Engel et al. (2013) points out that short periods of non-stationary network connectivity (termed intrinsic coupling modes; or ICMs) may occur between cortical regions either in the oscillatory signal or in the bandlimited envelope, and they are both potentially present in the task positive state. Further work is clearly needed to understand the mechanisms responsible for the phenomena observed here, and that should include measures of phase coherence between regions found to have high effective connectivity in the post-stimulus envelope. However, it is possible that measuring the post-stimulus envelope has allowed the measurement of some aspect of brain function that is not normally visible in the oscillatory phase -perhaps due to SNR issues in the oscillatory signal, or because those cortical processes that underlie these envelope modulations do not necessarily rely on phase coupling. We observed clear travelling wave like phenomena in the envelope that shares some similarities with TW phenomena in the oscillatory signal of the EEG and MEG.
It has been shown that cortical travelling exist in the oscillatory signal at varying spatial scales ranging from the columnar scale covering small patches of neurons a few millimetres across (Brookes et al., 2011a;Sato et al., 2012), intra-regional scales, such as within the entire primary visual cortex or motor cortex (Freeman and Barrie, 2000;Rubino et al., 2006;Takahashi et al., 2011) and inter-regional scale covering the whole brain (Alexander et al., 2006;Brookes et al., 2011aBrookes et al., , 2011aNunez, 1974a;Nunez et al., 1997). Nunez (1974a) introduced the concept of the global theory of standing and travelling waves and proposed the brain wave equation as a theoretical model for the global wave like phenomena observed in the EEG. The general postulate put forward (Nunez and Srinivasan, 2006) is that oscillatory waves observed in the EEG are partly composed of TWs, as the term is used in the physical sciences to describe a disturbance that moves through a medium (for example, as sound waves propagate through air). A number of experimental studies have demonstrated the existence of wave like behaviour in the EEG during wakefulness in the alpha band (Nunez, 1974b), slow wave sleep (Massimini et al., 2004(Massimini et al., , 2007(Massimini et al., , 2009, and in steady state visually evoked potentials (SSVEPs) (Burkitt et al., 2000;Hughes, 1995;Klimesch et al., 2007;Thorpe et al., 2007), demonstrating the general feasibility of viewing the cortex as a medium through which waves of excitatory and inhibitory action may travel. One major drawback in the experimental study of TWs is that they are degraded by trial averaging of the post stimulus response. In a study employing EEG, ECoG, and MEG, Alexander et al. (2011a) found that travelling waves analysed on a single trial basis accounted for the spatiotemporal pattern of observed ERs in trial-averaged data, but as evoked responses depend mainly on phase locking, trial averaging obscures travelling wave like phenomena. Evoked responses can therefore be explained as TWs and standing waves in the global field oscillations that are phase locked to the stimulus only at specific locations across the cortex, leading to an interference pattern in the trial averaged signal.
Consequently, while evoked responses may be a useful indicator of primary sensory processing latency (Henson and Rugg, 2003), it is likely that trial averaging obscures important information relating to cortical dynamics and stimulus driven connectivity. Those methods described above attempt to overcome this issue by measuring phase differences between adjacent sensors, building up a picture of the phase gradient (2011a, 2011a). The phase gradient therefore shows the general direction of wave propagation by indexing the phase coherence between adjacent sensors. However, phase gradients do not necessarily lead to long-range coherence as cumulative phase jitter may result in phase randomisation at long distances. Alexander found that the theta band TW moves in the posterioranterior direction while alpha moves in the anterior-posterior direction, which is broadly in agreement with both Bastos (Bastos et al., 2015a) and our results.
Research into travelling waves outlined above points to the possibility that the evoked response and the induced response are closely related; the induced response is simply a time shifted version of the evoked response with randomised phase in the carrier signal. In fact, there is some evidence to suggest that components of the evoked response that are lost in trial averaging of the raw signal due to phase variations may be transferred to the induced response (David et al., 2005;Makeig, 1993). These two aspects of post stimulus activity (evoked and induced) are therefore closely related and may together form a travelling wave like pattern in the envelope signal (that contains both evoked and induced responses) as the later cortical events occur at more distal regions from the origin. In our case, the evoked response in the visual cortex marks the origin of a travelling wave that becomes increasingly phase incoherent as it moves further away from the origin. It is therefore concluded that temporal delay observed between cortical sites involves the gradual transition from a phase locked to a nonphase locked response, which is only visible in the envelope signal. This idea is directly testable by measuring how the evoked signal relates to the travelling wave in the envelope observed here and could be the focus of further investigations.