Source Analysis 3D Imaging
Module information  
Modules  BESA Research Standard or higher 
Version  6.1 or higher 
Contents
 1 Overview
 2 Multiple Source Beamformer (MSBF) in the Timefrequency Domain
 3 Dynamic Imaging of Coherent Sources (DICS)
 4 Multiple Source Beamformer (MSBF) in the Time Domain
 5 CLARA
 6 LAURA
 7 LORETA
 8 sLORETA
 9 swLORETA
 10 sSLOFO
 11 UserDefined Volume Image
 12 Regularization of distributed volume images
 13 Cortical LORETA
 14 Cortical CLARA
 15 Cortex Inflation
 16 Surface Minimum Norm Image
 17 Multiple Source Probe Scan (MSPS)
 18 Source Sensitivity
 19 SESAME
 20 Brain Atlas
Overview
BESA Research features a set of new functions that provide 3D images that are displayed superimposed to the individual subject's anatomy. This chapter introduces these different images and describe their properties and applications.
The 3D images can be divided into three categories:
Volume images:
 The Multiple Source Beamformer (MSBF) is a tool for imaging brain activity. It is applied in the timedomain or timefrequency domain. The beamformer technique in timefrequency domain can image not only evoked, but also induced activity, which is not visible in timedomain averages of the data.
 Dynamic Imaging of Coherent Sources (DICS) can find coherence between any two pairs of voxels in the brain or between an external source and brain voxels. DICS requires timefrequencytransformed data and can find coherence for evoked and induced activity.
The following imaging methods provide an image of brain activity based on a distributed multiple source model:
 CLARA is an iterative application of LORETA images, focusing the obtained 3D image in each iteration step.
 LAURA uses a spatial weighting function that has the form of a local autoregressive function.
 LORETA has the 3D Laplacian operator implemented as spatial weighting prior.
 sLORETA is an unweighted minimum norm that is standardized by the resolution matrix.
 swLORETA is equivalent to sLORETA, except for an additional depth weighting.
 SSLOFO is an iterative application of standardized minimum norm images with consecutive shrinkage of the source space.
 A Userdefined volume image allows to experiment with the different imaging techniques. It is possible to specify userdefined parameters for the family of distributed source images to create a new imaging technique.
 Bayesian source imaging: SESAME uses a semiautomated Bayesian approach to estimate the number of dipoles along with their parameters.
Surface image:
 The Surface Minimum Norm Image. If no individual MRI is available, the minimum norm image is displayed on a standard brain surface and computed for standard source locations. If available, an individual brain surface is used to construct the distributed source model and to image the brain activity.
 Cortical LORETA. Unlike classical LORETA, cortical LORETA is not computed in a 3D volume, but on the cortical surface.
 Cortical CLARA. Unlike classical CLARA, cortical CLARA is not computed in a 3D volume, but on the cortical surface.
Discrete model probing:
These images do not visualize source activity. Rather, they visualize properties of the currently applied discrete source model:
 The Multiple Source Probe Scan (MSPS) is a tool for the validation of a discrete multiple source model.
 The Source Sensitivity image displays the sensitivity of a selected source in the current discrete source model and is therefore data independent.
Multiple Source Beamformer (MSBF) in the Timefrequency Domain
Short mathematical introduction
The BESA beamformer is a modified version of the linearly constrained minimum variance vector beamformer in the timefrequency domain as described in Gross et al., "Dynamic imaging of coherent sources: Studying neural interactions in the human brain", PNAS 98, 694699, 2001. It allows to image evoked and induced oscillatory activity in a userdefined timefrequency range, where time is taken relative to a triggered event.
The computation is based on a transformation of each channel's single trial data from the time domain into the timefrequency domain. This transformation is performed by the BESA Research Source Coherence module and leads to the complex spectral density S_{i} (f,t), where i is the channel index and f and t denote frequency and time, respectively. Complex cross spectral density matrices C are computed for each trial:
[math]\mathrm{C}_{ij}\left( f,t \right) = \mathrm{S}_{i}\left( f,t \right) \cdot \mathrm{S}_{j}^{*}\left( f,t \right)[/math]
The output power P of the beamformer for a specific brain region at location r is then computed by the following equation:
[math]\mathrm{P}\left( r \right) = \operatorname{tr^{'}}\left\lbrack \mathrm{L}^{T}\left( r \right) \cdot \mathrm{C}_{r}^{1} \cdot \mathrm{L}\left( r \right) \right\rbrack^{1}[/math]
Here, C_{r}^{1} is the inverse of the SVDregularized average of C_{ij}(f,t) over trials and the timefrequency range of interest; L is the leadfield matrix of the model containing a regional source at target location r and, optionally, additional sources whose interference with the target source is to be minimized; tr'[] is the trace of the [3×3] (MEG:[2×2]) submatrix of the bracketed expression that corresponds to the source at target location r.
In BESA Research, the output power P(r) is normalized with the output power in a reference timefrequency interval P_{ref}(r). A value q ist defined as follows:
[math] \mathrm{q}\left( r \right) =
\begin{cases}
\sqrt{\frac{\mathrm{P}\left( r \right)}{\mathrm{P}_{\text{ref}}(r)}}  1 = \sqrt{\frac{\operatorname{tr^{'}}\left\lbrack \mathrm{L}^{T}\left( r \right) \cdot \mathrm{C}_{r}^{ 1} \cdot \mathrm{L}\left( r \right) \right\rbrack^{ 1}}{\operatorname{tr^{'}}\left\lbrack \mathrm{L}^{T}\left( r \right) \cdot \mathrm{C}_{\text{ref},r}^{ 1} \cdot \mathrm{L}\left( r \right) \right\rbrack^{ 1}}}  1, & \text{for }\mathrm{P}(r) \geq \mathrm{P}_{\text{ref}}(r) \\
1  \sqrt{\frac{\mathrm{P}_{\text{ref}}\left( r \right)}{\mathrm{P}\left( r \right)}} = 1  \sqrt{\frac{\operatorname{tr^{'}}\left\lbrack \mathrm{L}^{T}\left( r \right) \cdot \mathrm{C}_{\text{ref},r}^{ 1} \cdot \mathrm{L}\left( r \right) \right\rbrack^{ 1}}{\operatorname{tr^{'}}\left\lbrack \mathrm{L}^{T}\left( r \right) \cdot \mathrm{C}_{r}^{ 1} \cdot \mathrm{L}\left( r \right) \right\rbrack^{ 1}}}, & \text{for }\mathrm{P}(r) \lt \mathrm{P}_{\text{ref}}(r)
\end{cases} [/math]
P_{ref }can be computed either from the corresponding frequency range in the baseline of the same condition (i.e. the beamformer images eventrelated power increase or decrease) or from the corresponding timefrequency range in a control condition (i.e. the beamformer images differences between two conditions). The beamformer image is constructed from values q(r) computed for all locations on a grid specified in the General Settings tab. For MEG data, the innermost grid points within a sphere of approx. 12% of the head diameter are assigned interpolated rather than calculated values).
qvalues are shown in %, where where q[%] = q*100. Alternatively to the definition above, q can also be displayed in units of dB:
[math]\mathrm{q}\left\lbrack \text{dB} \right\rbrack = 10 \cdot \log_{10}\frac{\mathrm{P}\left( r \right)}{\mathrm{P}_{\text{ref}}\left( r \right)}[/math]
A beamformer operator is designed to pass signals from the brain region of interest r without attenuation, while minimizing interference from activity in all other brain regions. Traditional singlesource beamformers are known to mislocalize sources if several brain regions have highly correlated activity. Therefore, the BESA beamformer extends the traditional singlesource beamformer in order to implicitly suppress activity from possibly correlated brain regions. This is achieved by using a multiple source beamformer calculation that contains not only the leadfields of the source at the location of interest r, but also those of possibly interfering sources. As a default, BESA Research uses a bilateral beamformer, where specifically contributions from the homologue source in the opposite hemisphere are taken into account (the matrix L thus being of dimension N×6 for EEG and N×4 for MEG, respectively, where N is the number of sensors). This allows for imaging of highly correlated bilateral activity in the two hemispheres that commonly occurs during processing of external stimuli.
In addition, the beamformer computation can take into account possibly correlated sources at arbitrary locations that are specified in the current solution. This is achieved by adding their leadfield vectors to the matrix L in the equation above.
Applying the Beamformer
This chapter illustrates the usage of the BESA beamformer. The displayed figures are generated using the file 'Examples/LearnbySimulations/ACCoherence/ACOsc20.foc' (see BESA Tutorial 6: "Timefrequency analysis and Source coherence").
Starting the beamformer from the timefrequency window
The BESA beamformer is applied in the timefrequency domain and therefore requires the Source Coherence module to be enabled. The timefrequency beamformer is especially useful to image in or decrease of induced oscillatory activity. Induced activity cannot be observed in the averaged data, but shows up as enhanced averaged power in the TSE (TemporalSpectral Evolution) plot. For instructions on how to initiate a beamformer computation in the timefrequency window, please refer to Chapter How to Create Beamformer Images.
After the beamformer computation has been initiated in the timefrequency window, the source analysis window opens with an enlarged 3D image of the qvalue computed with a bilateral beamformer. The result is superimposed onto the MR image assigned to the data set (individual or standard).
Beamformer image after starting the computation in the TimeFrequency window. A bilateral pair of sources in the auditory cortex accounts for the highly correlated oscillatory induced activity. Only the bilateral beamformer manages to separate these activities; a traditional singlesource beamformer would merge the two sources into one image maximum in the head center instead.
Multiple source beamformer in the Source Analysis window
The 3D imaging display is part of the source analysis window. If you press the Restore button at the right end of the title bar of the 3D window, the window appears at the bottom right of the source analysis window. In the channel box, the averaged (evoked) data of the selected condition is shown. When a control condition was selected, its average is appended to the average of the target condition.
Source Analysis window with beamformer image. The two sources have been added using the Switch to Maximum and Add Source toolbar buttons (see below). Source waveforms are computed from the displayed averaged data. Therefore, they do not represent the activity displayed in the beamformer image, which in this simulation example is induced (i.e. not phaselocked to the trigger)!
When starting the beamformer from the timefrequency window, a bilateral beamformer scan is performed. In the source analysis window, the beamformer computation can be repeated taking into account possibly correlated sources that are specified in the current solution. Interfering activities generated by all sources in the current solution that are in the 'On' state are specifically suppressed (they enter the matrix L in the beamformer calculation, see Chapter Short mathematical description above). The computation can be started from the Image menu or from the Image selector button dropdown menu. The Image menu can be evoked either from the menu bar or by rightclicking anywhere in the source analysis window.
Multiple source beamformer image calculated in the presence of a source in the left hemisphere. A single source scan has been performed. The source set in the current solution accounts for the lefthemispheric qmaximum in the data. Accordingly, the beamformer scan reveals only the as yet unmodeled additional activity in the right hemisphere (note the radiological convention in the 3D image display).
The beamformer scan can be performed with a single or a bilateral source scan. The default scan type depends on the current solution:
 When the beamformer is started from the TimeFrequency window, the Source Analysis window opens with a new solution and a bilateral beamformer scan is performed.
 When the beamformer is started within the Source Analysis window, the default is
 a scan with a single source in addition to the sources in the current solution, if at least one source is active.
 a bilateral scan if no source in the current solution is active.
The default scan type is the multiple source beamformer. The nondefault scan type can be enforced using the corresponding Volume Image / Beamformer entry in the Image menu.
Inserting Sources out of the Beamformer Image
The beamformer image can be used to add sources to the current solution. A simple doubleclick anywhere in the 2D or 3Dview will generate a nonoriented regional source at the corresponding location. However, a better and easier way to create sources at image maxima and minima is to use the toolbar buttons Switch to Maximum and Add Source .
Use the Switch to Maximum button to place the red crosshair of the 3D window onto a local image maximum or minimum. Hitting the Add Source button creates a regional source at the location of the crosshair and therefore ensures the exact placement of the source at the image extremum. Moreover, the Add Source button generates an oriented regional source. BESA Research automatically estimates the source orientation that contributes most to the power in the target timefrequency interval (or the reference timefrequency interval, if its power is larger than that in the target interval). The accuracy of this orientation estimate depends largely on the noise content of the data. The smaller the signaltonoise ratio of the data, the lower is the accuracy of the orientation estimate. This feature allows to use the beamformer as a tool to create a source montage for source coherence analysis, where it is of advantage to work with oriented sources.
Notes:
 You can hide or redisplay the last computed image by selecting the corresponding entry in the Image menu.
 The current image can be exported to ASCII or BrainVoyager vmpformat from the Image menu.
 For scaling options, use the and Image Scale toolbar buttons.
 Parameters used for the beamformer calculations can be set in the Standard Volumes of the Image Settings dialog box.
Dynamic Imaging of Coherent Sources (DICS)
Short mathematical introduction
Dynamic Imaging of Coherent Sources (DICS) is a sophisticated method for imaging corticocortical coherence in the brain, or coherence between an external reference (e.g. EMG channel) and cortical structures. DICS can be applied to localize evoked as well as induced coherent cortical activity in a userdefined timefrequency range.
DICS was implemented in BESA closely following Gross et al., "Dynamic imaging of coherent sources: Studying neural interactions in the human brain", PNAS 98, 694699, 2001.
The computation is based on a transformation of each channel's single trial data from the time domain into the frequency domain. This transformation is performed by the BESA Research Coherence module and results in the complex spectral density matrix that is used for constructing the spatial filter similar to beamforming.
DICS computation yields a 3D image, each voxel being assigned a coherence value. Coherence values can be described as a neural activity index and do not have a unit. The neural activity index contrasts coherence in a target timefrequency bin with coherence of the same timefrequency bin in a baseline.
DICS for corticocortical coherence is computed as follows:
Let L(r) be the leadfield in voxel r in the brain and C the complex crossspectral density matrix. The spatial filter W(r) for the voxel r in the head is defined as follows:
[math]W\left( r \right) = \left\lbrack L^{T}\left( r \right) \cdot C^{ 1} \cdot L\left( r \right) \right\rbrack^{ 1} \cdot L^{T}(r) \cdot C^{ 1}[/math]
The crossspectrum between two locations (voxels) r_{1} and r_{2} in the head are calculated with the following equation:
[math]C_{s}\left( r_{1},r_{2} \right) = W\left( r_{1} \right) \cdot C \cdot W^{*T}\left( r_{2} \right),[/math]
where ^{*T} means the transposed complex conjugate of a matrix. The crossspectral density can then be calculated from the cross spectrum as follows:
[math]c_{s}\left( r_{1},r_{2} \right) = \lambda_{1}\left\{ C_{s}\left( r_{1},r_{2} \right) \right\},[/math]
where λ_{1}{} indicates the largest singular value of the cross spectrum. Once the cross spectral density is estimated, the connectivity¹(CON) between the two brain regions r_{1} and r_{2} are calculated as follows:
[math]\text{CON}\left( r_{1},r_{2} \right) = \frac{c_{s}^{\text{sig}}\left( r_{1},r_{2} \right)  c_{s}^{\text{bl}}(r_{1},r_{2})}{c_{s}^{\text{sig}}\left( r_{1},r_{2} \right) + c_{s}^{\text{bl}}(r_{1},r_{2})},[/math]
where c_{s}^{sig} is the crossspectral density for the signal of interest between the two brain regions r_{1} and r_{2}, and c_{s}^{bl} is the corresponding cross spectral density for the baseline or the control condition, respectively. In the case DICS is computed with a cortical reference, r_{1} is the reference region (voxel) and remains constant while r_{2} scans all the grid points within the brain sequentially. In that way, the connectivity between the reference brain region and all other brain regions is estimated. The value of CON(r_{1}, r_{2}) falls in the interval [1 1]. If the crossspectral density for the baseline is 0 the connectivity value will be 1. If the crossspectral density for the signal is 0 the connectivity value will be 1.
¹ Here, the term connectivity is used rather than coherence, as strictly speaking the coherence equation is defined slightly differently. For simplicity reasons the rest of the tutorial uses the term coherence.
DICS for corticomuscular coherence is computed as follows:
When using an external reference, the equation for coherence calculation is slightly different compared to the equation for corticocortical coherence. First of all, the crossspectral density matrix is not only computed for the MEG/EEG channels, but the external reference channel is added. This resulting matrix is C_{all}. In this case, the crossspectral density between the reference signal and all other MEG/EEG
channels is called c_{ref}. It is only one column of C_{all}. Hence, the crossspectrum in voxel r is calculated with the following equation:
[math]C_{s}\left( r \right) = W\left( r \right) \cdot c_{\text{ref}}[/math]
and the corresponding crossspectral density is calculated as the sum of squares of C_{s}:
[math]c_{s}\left( r \right) = \sum_{i = 1}^{n}{C_{s}\left( r \right)_{i}^{2}},[/math]
where n is 2 for MEG and 3 for EEG. This equation can also be described as the squared Euclidean norm of the crossspectrum:
[math]c_{s}\left( r \right) = \left\ C_{s} \right\^{2},[/math]
The power in voxel r is calculated as in the corticocortical case:
[math]p\left( r \right) = \lambda_{1}\left\{ C_{s}(r,r) \right\}.[/math]
At last, coherence between the external reference and cortical activity is calculated with the equation:
[math]\text{CON}\left( r \right) = \frac{c_{s}(r)}{p\left( r \right) \cdot C_{\text{all}}(k,k)},[/math]
where C_{all}(k, k) is the (k,k)th diagonal element of the matrix C_{all}.
DICS is particularly useful, if coherence is to be calculated without an apriory source model (in contrast to source coherence based on predefined source montages). However, the recommended analysis strategy for DICS is to use a brain source as a starting point for coherence calculation that is known to contribute to the EEG/MEG signal of interest. For example, one might first run a beamformer on the timefrequency range of interest and use the voxel with the strongest oscillatory activity as a starting point for DICS. The resulting coherence image will again lead to several maxima (ordered by magnitude), which in turn can serve as starting points for DICS calculation. This way, it is possible to detect even weak sources that show coherent activity in the given timefrequency range.
The other significant application for DICS is estimating coherence between an external source and voxels in the brain. For example, an external source can be muscle activity recoded by an electrode placed over the according peripheral region. This way, the direct relationship between muscle activity and brain activation can be measured.
Starting DICS computation from the TimeFrequency Window
DICS is particularly useful, if coherence in a userdefined timefrequency bin (evoked or induced) is to be calculated between any two brain regions or between an external reference and the brain. DICS runs only on timefrequency decomposed data, so timefrequency analysis needs to be run before starting DICS computation.
To start the DICS computation, leftdrag a window over a selected timefrequency bin in the TimeFrequency Window. Rightclick and select “Image”. A dialogue will open (see fig. 1) prompting you to specify time and frequency settings as well as the baseline period. It is recommended to use a baseline period of equal length as the data period of interest. Make sure to select “DICS” in the top row and press “Go”.
Next, a window will appear allowing you to specify the reference source for coherence calculation (see fig. 2). It is possible to select a channel (e.g. EMG) or a brain source. If a brain source is chosen and no source analysis was computed beforehand, the option “Use current crosshair position” must be chosen. In case discrete source analysis was computed previously, the selected source can be chosen as the reference for DICS. Please note that DICS can be recomputed with any crosshair or source position at a later stage.
Confirming with “OK” will start computation of coherence between the selected channel/voxel and all other brain voxels. In case DICS is computed for a reference source in the brain, it can be advantageous to run a beamforming analysis in the selected timefrequency window first and use one of the beamforming maxima as reference for DICS. Fig. 3 shows an example for DICS calculation.
Coherence values range between 1 and 1. If coherence in the signal is much larger than coherence in the baseline (control condition) then the DICS value is going to approach 1. Contrary, if coherence in the baseline is much larger than coherence in the signal, then the DICS value is going to approach 1. At last, if coherence in the signal is equal to coherence in the baseline, then the DICS value is 0.
In case DICS is to be recomputed with a different reference, simply mark the desired reference position by placing the crosshair in the anatomical view and select “DICS” in the middle panel of the source analysis window (see Fig. 4). In case an external reference is to be selected, click on “DICS” in the middle panel to bring up the DICS dialogue (see. Fig. 2) and select the desired channel. Please note that DICS computation will only be available after running timefrequency analysis.
Multiple Source Beamformer (MSBF) in the Time Domain
(requires Besa Research 7.0 or higher)
Short mathematical introduction
Beamforming approach can be also applied in the time domain data. This approach was introduced as linearly constrained minimum variance (LCMV) beamformer (Van Veen et al., 1997). It allows to image evoked activity in a userdefined time range, where time is taken relative to a triggered event, and to estimate source waveforms using the calculated spatial weight at locations of interest. For an implementation of the beamformer in the time domain, data covariance matrices are required, while complex cross spectral density matrices are used for the beamformer approaches in the timefrequency domain as described in the Multiple Source Beamformer (MSBF) in the Timefrequency Domain section.
The bilateral beamformer introduced in the Multiple Source Beamformer (MSBF) in the Timefrequency Domain section is also implemented for the timedomain beamformer to take into account contributions from the homologue source in the opposite. This allows for imaging of highly correlated bilateral activity in the two hemispheres that commonly occurs during processing of external stimuli. In addition, the beamformer computation can take into account possibly correlated sources at arbitrary locations.
The beamformer spatial weight W(r) for the voxel r in the brain is defined as follows (Van Veen et al., 1997):
[math]\mathrm{W}(r) = [\mathrm{L}^T(r)\mathrm{C}^{1}\mathrm{L}(r)]^{1}\mathrm{L}^T(r)\mathrm{C}^{1}[/math]
where [math]\mathrm{C}^{1}[/math] is the inversed regularized average of covariance matrix over trials, L is the leadfield matrix of the model containing a regional source at target location r and optionally additional sources whose interference with the target source is to be minimized. The beamformer spatial weight W(r) can be applied to the measured data to estimate source waveform at a location r (beamformer virtual sensor):
[math]\mathrm{S}(r,t) = \mathrm{W}(r)\mathrm{M}(t)[/math]
where S(r,t) represents the estimated source waveform and M(t) represents measured EEG or MEG signals. The output power P of the beamformer for a specific brain region at location r is computed by the following equation:
[math]\mathrm{P}(r) = \operatorname{tr^{'}}[\mathrm{W}(r) \cdot \mathrm{C} \cdot \mathrm{W}^T(r)][/math]
where tr’[ ] is the trace of the [3×3] (MEG: [2×2]) submatrix of the bracketed expression that corresponds to the source at target location r.
Beamformer can suppress noise sources that are correlated across sensors. However, uncorrelated noise will be amplified in a spatially nonuniform manner, with increasing distortion with increasing distance from the sensors (Van Veen et al., 1997; Sekihara et al., 2001). For this reason, estimated source power should be normalized by a noise power. In BESA Research, the output power P(r) is normalized with the output power in a baseline interval or with the output power of a uncorrelated noise: P(r) / Pref (r).
The timedomain beamformer image is constructed from values q(r) computed for all locations on a grid specified in the General Settings tab. A value q(r) is defined as described in the Multiple Source Beamformer (MSBF) in the Timefrequency Domain section with data covariance matrices instead of crossspectral density matrices.
Applying the Beamformer
This chapter illustrates the usage of the BESA beamformer in the time domain. The displayed figures are generated using the file ‘Examples/ERPAuditoryIntensity/S1.cnt’.
Starting the timedomain beamformer from the Average tab of the Paradigm dialog box
The timedomain beamformer is needed data covariance matrices and therefore requires the ERP module to be enabled. After the beamformer computation has been initiated in the Average tab of the Paradigm dialog box, the source analysis window opens with an enlarged 3D image of the qvalue computed with a bilateral beamformer. The result is superimposed onto the MR image assigned to the data set (individual or standard).
Multiplesource beamformer in the Source Analysis window
The 3D imaging display is part of the source analysis window. In the Channel box, the averaged (evoked) data of the selected condition is shown. Selected covariance intervals in the ERP module can be checked in the Channel box. The red, gray, and blue rectangles indicate signal, baseline, and common interval, respectively.
When starting the beamformer from the Average tab of the Paradigm dialog box, the bilateral beamformer scan is performed. In the source analysis window, the beamformer computation can be repeated taking into account possibly correlated sources that are specified in the current solution. Interfering activities generated by all sources in the current solution that are in the 'On' state are specifically suppressed (they enter the leadfield matrix L in the beamformer calculation). The computation can be started from the Image menu or from the Image selector button dropdown menu. The Image menu can be evoked either from the menu bar or by rightclicking anywhere in the source analysis window.
The beamformer scan can be performed with a single or a bilateral source scan. The default scan type depends on the current solution:
When the beamformer is started from the Average tab of the Paradigm dialog box the Source Analysis window opens with a new solution and a bilateral beamformer scan is performed.
When the beamformer is started within the Source Analysis window, the default is:
 a scan with a single source in addition to the sources in the current solution, if at least one source is active.
 a bilateral scan if no source in the current solution is active.
 a scan with a single source when scalartype beamformer is selected in the beamformer option dialog box.
The default scan type is the multiple source beamformer. The nondefault scan type can be enforced using the corresponding Volume Image / Beamformer entry in the Image main menu or in the beamformer option dialog box (only for the timedomain beamformer).
Inserting Sources as Beamformer Virtual Sensor out of the Beamformer Image
This is similar to the inserting sources out of the beamformer image in Multiple Source Beamformer (MSBF) in the Timefrequency Domain section.
The beamformer image can be used to add beamformer virtual sensors to the current solution. A simple doubleclick anywhere in the 3D view (not in the 2D view) will generate a source at the corresponding location. A better and easier way to create sources at image maxima and minima is to use the toolbar buttons Switch to Maximum and Add Source .
This feature allows to use the beamformer as a tool to create a source montage for source coherence analysis. A source montage file (*.mtg) for beamformer virtual sensors can be saved using File \ Save Source Montage As… entry.
The timedomain beamformer image can be also used to add regional or dipole sources to the current solution. Press N key when there is no source in the current source array or there is more than one beamformer virtual sensor. To create a new source array for beamformer virtual sensor, press N key when there is more than one regional or dipole source in the current source array.
Notes
 You can hide or redisplay the last computed image by selecting Hide Image entry in the Image menu.
 The current image can be exported to ASCII, ANALYZE, or BrainVoyager (*.vmp) format from the Image menu.
 For scaling options, use and Image Scale toolbar buttons.
 Parameters used for the beamformer calculations can be set in the Standard Volume tab of the Image Settings dialog box.
 Note that Model, Residual, Order, and Residual variance are not shown for the beamformer virtual sensor type sources.
References
 Sekihara, K., Nagarajan, S. S., Poeppel, D., Marantz, A., & Miyashita, Y. (2001). Reconstructing spatiotemporal activities of neural sources using an MEG vector beamformer technique. IEEE Transactions on Biomedical Engineering, 48(7), 760–771.
 Van Veen, B. D., Van Drongelen, W., Yuchtman, M., & Suzuki, A. (1997). Localization of brain electrical activity via linearly constrained minimum variance spatial filtering. IEEE Transactions on Biomedical Engineering, 44(9), 867–880
CLARA
CLARA ('Classical LORETA Analysis Recursively Applied') is an iterative application of weighted LORETA images with a reduced source space in each iteration.
In an initialization step, a LORETA image is calculated. Then in each iteration the following steps are performed:
 The obtained image is spatially smoothed (this step is left out in the first iteration).
 All grid points with amplitudes below a threshold of 1% of the maximum activity are set to zero, thus being effectively eliminated from the source space in the following step.
 The resulting image defines a spatial weighting term (for each voxel the corresponding image amplitude).
 A LORETA image is computed with an additional spatial weighting term for each voxel as computed in step 3. By the default settings in BESA Research, the regularization values used in the iteration steps are slightly higher than that of the initialization LORETA image.
The procedure stops after 2 iterations, and the image computed in the last iteration is displayed. Please note that you can change all parameters by creating a userdefined volume image.
The advantage of CLARA over nonfocusing distributed imaging methods is visualized by the figure below. Both images are computed from the N100 response in an auditory oddball experiment (file Oddball.fsg in subfolder fMRI+EEGRTExperiment of the Examples folder). The CLARA image is much more focal than the sLORETA image, making it easier to determine the location of the image maxima.
Notes:
 Starting CLARA: CLARA can be started from the Image menu or from the Image Selection button.
 Please refer to Chapter Regularization of distributed volume images for important information on regularization of distributed inverses.
LAURA
LAURA (Local Auto Regressive Average) belongs to the distributed inverse method of the family of weighted minimum norm methods (Grave de Peralta Menendeza et al., "Noninvasive Localization of Electromagnetic Epileptic Activity. I. Method Descriptions and Simulations", BrainTopography 14(2), 131137, 2001). LAURA uses a spatial weighting function that includes depth weighting and that term has the form of a local autoregressive function.
The source activity is estimated by applying the general formula for a weighted minimum norm:
[math]\mathrm{S}\left( t \right) = \mathrm{V} \cdot \mathrm{L}^{T}\left( \mathrm{L} \cdot \mathrm{V} \cdot \mathrm{L}^{T} \right)^{ 1} \cdot \mathrm{D}(t)[/math]
Here, L is the leadfield matrix of the distributed source model with regional sources distributed on a regular cubic grid. D(t) is the data at time point t. The term in parentheses is generally regularized. Regularization parameters can be specified in the Image Settings.
In LAURA, V contains both a depth weighting term W and a representation of a local autoregressive function A. V is computed as:
[math]\mathrm{V} = \left( \mathrm{U}^{T} \cdot \mathrm{U} \right)^{1}[/math]
Where
[math]\mathrm{U} = \left( \mathrm{W} \cdot \mathrm{A} \right) \otimes \mathrm{I}_{3}[/math]
Here, [math]\otimes[/math] denotes the Kronecker product. I_{3} is the [3×3] identity matrix. W is an [s×s] diagonal matrix (with s the number of source locations on the grid), where each diagonal element is the inverse of the maximum singular value of the corresponding regional source's leadfields. The formula for the diagonal components A_{ii} and the offdiagonal components A_{ik} are as follows:
[math]\mathrm{A}_{ii} = \frac{26}{\mathrm{N}_{i}}\sum_{k \subset V_{i}}^{}\frac{1}{\mathrm{d}_{ik}^{2}}[/math]
[math] \mathrm{A}_{ik} = \begin{cases}  1/\operatorname{dist}\left( i,k \right)^{2}, & \text{if } k \subset V_{i} \\ 0, & \text{otherwise} \end{cases} [/math]
Here, V_{i} is the vicinity around grid point i that includes the 26 direct neighbors.
The LAURA image in BESA Research displays the norm of the 3 components of S at each location r. Using the menu function Image / Export Image As... you have the option to save this norm of S or alternatively all components separately to disk.
Notes:
 Grid spacing: Due to memory limitations, LAURA images require a grid spacing of 7 mm or more.
 Computation time: Computation speed during the first LAURA image calculation depends on the grid spacing (computation is faster with larger grid spacing). After the first computation of a LAURA image, a *.laura file is stored in the data folder, containing intermediate results of the LAURA inverse. This file is used during all subsequent LAURA image computations. Thereby, the time needed to obtain the image is substantially reduced.
 MEG: In the case of MEG data, an additional constraint is implemented in the LAURA algorithm that prevents solutions from containing radial source currents (compare PascualMarqui, ISBET Newsletter 1995, 2229). In MEG, an additional source space regularization is necessary in the inverse matrix operation required compute V
 Starting LAURA: LAURA can be started from the Image menu or from the Image Selection button.
 Regularization: Please refer to Chapter “Regularization of distributed volume images” for important information on regularization of distributed inverses.
LORETA
LORETA ("Low Resolution Electromagnetic Tomography") is a distributed inverse method of the family of weighted minimum norm methods. LORETA was suggested by R.D. PascualMarqui (International Journal of Psychophysiology. 1994, 18:4965). LORETA is characterized by a smoothness constraint, represented by a discrete 3D Laplacian.
The source activity is estimated by applying the general formula for a weighted minimum norm:
[math]\mathrm{S}\left( t \right) = \mathrm{V} \cdot \mathrm{L}^{T}\left( \mathrm{L} \cdot \mathrm{V} \cdot \mathrm{L}^{T} \right)^{ 1} \cdot \mathrm{D}(t)[/math]
Here, L is the leadfield matrix of the distributed source model with regional sources distributed on a regular cubic grid. D(t) is the data at time point t. The term in parentheses is generally regularized. Regularization parameters can be specified in the Image Settings.
In LORETA, V contains both a depth weighting term and a representation of the 3D Laplacian matrix. V is computed as:
[math]\mathrm{V} = \left( \mathrm{U}^{T} \cdot \mathrm{U} \right)^{ 1}[/math]
where
[math]\mathrm{U} = \left( \mathrm{W} \cdot \mathrm{A} \right) \otimes \mathrm{I}_{3}[/math]
Here, [math]\otimes[/math] denotes the Kronecker product. I_{3} is the [3x3] identity matrix. W is an [sxs] diagonal matrix (with s the number of source locations on the grid), where each diagonal element is the inverse of the maximum singular value of the corresponding regional source's leadfields. A contains the 3D Laplacian and is computed as
[math]\mathrm{A} = \mathrm{Y}  \mathrm{I}_{s}[/math]
with I_{s} the [sxs] identity matrix, where s is the number of sources (= three times the number of grid points) and
[math]\mathrm{Y} = \frac{1}{2}\left\{ \mathrm{I}_{s} + \left\lbrack \operatorname{diag}\left( \mathrm{Z} \cdot \left\lbrack 111 \ldots 1 \right\rbrack^{T} \right) \right\rbrack^{ 1} \right\} \cdot \mathrm{Z}[/math]
where
[math]
\mathrm{Z}_{ik} =
\begin{cases}
1/6, & \text{if } \operatorname{dist}\left( i,k \right) = 1 \text{ grid point} \\
0, & \text{otherwise}
\end{cases}
[/math]
The LORETA image in BESA Research displays the norm of the 3 components of S at each location r. Using the menu function Image / Export Image As... you have the option to save this norm of S or alternatively all components separately to disk.
Notes:
 Grid spacing: Due to memory limitations, LORETA images require a grid spacing of 5 mm or more.
 Computation time: Computation speed during the first LORETA image calculation depends on the grid spacing (computation is faster with larger grid spacing). After the first computation of a LORETA image, a *.loreta file is stored in the data folder, containing intermediate results of the LORETA inverse. This file is used during all subsequent LORETA image computations. Thereby, the time needed to obtain the image is substantially reduced.
 MEG: In the case of MEG data, an additional constraint is implemented in the LORETA algorithm that prevents solutions from containing radial source currents (PascualMarqui, ISBET Newsletter 1995, 2229). In MEG, an additional source space regularization is necessary in the inverse matrix operation required compute V.
 Starting LORETA: LORETA can be started from the Image menu or from the Image Selection button.
 Regularization: Please refer to Chapter “Regularization of distributed volume images” for important information on regularization of distributed source models.
sLORETA
This distributed inverse method consists of a standardized, unweighted minimum norm. The method was originally suggested by R.D. PascualMarqui (Methods & Findings in Experimental & Clinical Pharmacology 2002, 24D:512) Starting point is an unweighted minimum norm computation:
[math]\mathrm{S}_{\text{MN}}\left( t \right) = \mathrm{L}^{T} \cdot \left( \mathrm{L} \cdot \mathrm{L}^{T} \right)^{ 1} \cdot \mathrm{D}(t)[/math]
Here, L is the leadfield matrix of the distributed source model with regional sources distributed on a regular cubic grid. D(t) is the data at time point t. The term in parentheses is generally regularized. Regularization parameters can be specified in the Image Settings.
This minimum norm estimate is now standardized to produce the sLORETA activity at a certain brain location r:
[math]\mathrm{S}_{\text{sLORETA}, r} = \mathrm{R}_{rr}^{1/2} \cdot \mathrm{S}_{\text{MN},r}[/math]
S_{sMN,r }is the [3x1] (MEG: [2x1]) minimum norm estimate of the 3 (MEG: 2) dipoles at location r. R_{rr} is the [3x3] (MEG: [2x2]) diagonal block of the resolution matrix R that corresponds to the source components at the target location r. The resolution matrix is defined as:
[math]\mathrm{R} = \mathrm{L}^{T} \cdot \left( \mathrm{L} \cdot \mathrm{L}^{T} + \lambda \cdot \mathrm{I} \right)^{1} \cdot \mathrm{L}[/math]
The sLORETA image in BESA Research displays the norm of S_{sLORETA}, _{r} at each location r. Using the menu function Image / Export Image As... you have the option to save this norm of S_{sLORETA}, _{r} or alternatively all components separately to disk.
Notes:
 sLORETA can be started from the Image menu or from the Image Selection button.
 Please refer to Chapter Regularization of distributed volume images for important information on regularization of distributed inverses.
swLORETA
This distributed inverse method is a standardized, depthweighted minimum norm (E. PalmeroSoler et al 2007 Phys. Med. Biol. 52 17831800). It differs from sLORETA only by an additional depth weighting.
Starting point is a depthweighted minimum norm computation:
[math]\mathrm{S}_{\text{MN}}\left( t \right) = \mathrm{V} \cdot \mathrm{L}^{T} \cdot \left( \mathrm{L} \cdot \mathrm{V} \cdot \mathrm{L}^{T} \right)^{1} \cdot \mathrm{D}(t)[/math]
Here, L is the leadfield matrix of the distributed source model with regional sources distributed on a regular cubic grid. D(t) is the data at time point t. The term in parentheses is generally regularized. Regularization parameters can be specified in the Image Settings.
V is the diagonal depth weighting matrix. For s grid locations, V is of dimension [3s x 3s] (MEG: [2s x 2s]). Each diagonal element of V is the inverse of the first singular value of the leadfield of the corresponding regional source. Hence, the first 3 (MEG: 2) diagonal elements equal the inverse of the largest eigenvalue of the leadfield matrix of regional source 1, and so on.
This minimum norm estimate is now standardized to produce the swLORETA activity at a certain brain location r:
[math]\mathrm{S}_{\text{swLORETA},r} = \mathrm{R}_{rr}^{1/2} \cdot \mathrm{S}_{\text{MN},r}[/math]
S_{sMN,r} is the [3x1] (MEG: [2x1]) depthweighted minimum norm estimate of the regional source at location r. R_{rr} is the [3x3] (MEG: [2x2]) diagonal block of the resolution matrix R that corresponds to the source components at the target location r. The resolution matrix is defined as:
[math]\mathrm{R} = \mathrm{V} \cdot \mathrm{L}^{T} \cdot \left( \mathrm{L} \cdot \mathrm{V} \cdot \mathrm{L}^{T} + \lambda \cdot \mathrm{I} \right)^{1} \cdot \mathrm{L}[/math]
The swLORETA image in BESA Research displays the norm of S_{swLORETA}, _{r} at each location r. Using the menu function Image / Export Image As... you have the option to save this norm of S_{swLORETA}, _{r} or alternatively all components separately to disk.
Notes:
 sLORETA can be started from the Image menu or from the Image Selection button.
 Please refer to Chapter “Regularization of distributed volume images” for important information on regularization of distributed inverses.
sSLOFO
SSLOFO (standardized shrinking LORETAFOCUSS) is an iterative application of weighted distributed source images with a reduced source space in each iteration (Liu et al., "Standardized shrinking LORETAFOCUSS (SSLOFO): a new algorithm for spatiotemporal EEG source reconstruction", IEEE Transactions on Biomedical Engineering 52(10), 16811691, 2005).
In an initialization step, an sLORETA image is calculated. Then in each iteration the following steps are performed:
 A weighted minimum norm solution is computed according to the formula [math]\mathrm{S} = \mathrm{V} \cdot \mathrm{L}^{T} \cdot \left( \mathrm{L} \cdot \mathrm{V} \cdot \mathrm{L}^{T} \right)^{1} \cdot \mathrm{D}[/math] . Here, L is the leadfield matrix of the distributed source model with regional sources distributed on a regular cubic grid. D is the data at the time point under consideration. V is a diagonal spatial weighting matrix that is computed in the previous iteration step. In the first iteration, the elements of V contain the magnitudes of the initially computed LORETA image.
 Standardization of this weighted minimum norm image is performed with the resolution matrix as in sLORETA.
 The obtained standardized weighted minimum norm image is being smoothed to get S_{smooth}.
 All voxels with amplitudes below a threshold of 1% of the maximum activity get a weight of zero in the next iteration step, thus being effectively eliminated from the source space in the next iteration step.
 For all other voxels, compute the elements of the spatial weighting matrix V to be used in the next iteration as follows: [math]\mathrm{V}_{ii,\text{next iteration}} = \frac{1}{\left\ \mathrm{L}_{i} \right\} \cdot \mathrm{S}_{ii,\text{smooth}} \cdot \mathrm{V}_{ii,\text{current iteration}}[/math]
The procedure stops after 3 iterations. Please note that you can change all parameters by creating a userdefined volume image.
Notes:
 Starting sSLOFO: sSLOFO can be started from the Image menu or from the Image Selection button.
 Please refer to Chapter Regularization of distributed volume images for important information on regularization of distributed inverses.
UserDefined Volume Image
In addition to the predefined 3D imaging methods in BESA Research, it is possible to create userdefined imaging methods based on the general formula for distributed inverses:
[math]\mathrm{S}\left( t \right) = \mathrm{V} \cdot \mathrm{L}^{T} \cdot \left( \mathrm{L} \cdot \mathrm{V} \cdot \mathrm{L}^{T} \right)^{1} \cdot \mathrm{D}(t)[/math]
Here, L is the leadfield matrix of the distributed source model with regional sources distributed on a regular cubic grid. D(t) is the data at time point t. Customdefined parameters are:
 The spatial weighting matrix V: This may include depth weighting, image weighting, or crossvoxel weighting with a 3D Laplacian (as in LORETA) or an autoregressive function (as in LAURA).
 Regularization: The term in parentheses is generally regularized. Note that regularization has a strong effect on the obtained results. Please refer to chapter Regularization of Distributed Volume Imagesfor more information.
 Standardization: Optionally, the result of the distributed inverse can be standardized with the resolution matrix (as in sLORETA).
 Iterations: Inverse computations can be applied iteratively. Each iteration is weighted with the image obtained in the previous iteration.
All parameters for the userdefined volume image are specified in the UserDefined Volume Tab of the Image Settings dialog box. Please refer to chapter UserDefined Volume Tab for details.
Notes:
 Starting the userdefined volume image: the image calculation can be started from the Image menu or from the Image Selection button.
 Please refer to Chapter Regularization of distributed volume images for important information on regularization of distributed inverses.
Regularization of distributed volume images
Distributed source images require the inversion of a term of the form L V L^{T}. This term is generally regularized before its inversion. In BESA Research, selection can be made between two different regularization approaches (parameters are defined in the Image Settings dialog box):
 Tikhonov regularization: In Tikhonov regularization, the term L V L^{T} is inverted as (L V L^{T }+λ I)^{1}. Here, l is the regularization constant, and I is the identity matrix.
 One way of determining the optimum regularization constant is by minimizing the generalized cross validation error (CVE).
 Alternatively, the regularization constant can be specified manually as a percentage of the trace of the matrix L V L^{T}.
 TSVD: In the truncated singular value decomposition (TSVD) approach, an SVD decomposition of L V L^{T} is computed as L V L^{T} = U S U^{T}, where the diagonal matrix S contains the singular values. All singular values smaller than the specified percentage of the maximum singular values are set to zero. The inverse is computed as U S^{1} U^{T}, where the diagonal elements of S^{1 }are the inverse of the corresponding nonzero diagonal elements of S.
Regularization has a critical effect on the obtained distributed source images. The results may differ completely with different choices of the regularization parameter (see examples below). Therefore, it is important to evaluate the generated image critically with respect to the regularization constant, and to keep in mind the uncertainties resulting from this fact when interpreting the results. The default setting in BESA Research is a TSVD regularization with a 0.03% threshold. However, this value might need to be adjusted to the specific data set at hand.
The following example illustrates the influence of the regularization parameter on the obtained images. The data used here is condition StCor of dataset Examples \ TFCErrorRelatedNegativity \ Correct+Error.fsg at 176 ms following the visual stimulus. Discrete dipole analysis reveals the main activity in the left and right lateral visual cortex at this latency.
LORETA images computed at this latency depend critically on the choice of the regularization constant. The following 3D images are created with TSVD regularization with SVD cutoffs of 0.1%, 0.005%, and 0.0001%, respectively. The volume grid size was 9 mm. The example demonstrates the dramatic effect of regularization and demonstrates the typical tradeoff between too strong regularization (leading to too smeared 3D images that tend to show blurred maxima) and too small regularization (resulting in too superficial 3D images with multiple maxima).
The automatic determination of the regularization constant using the CVE approach does not necessarily result in the optimum regularization parameter either. In this example, the unscaled CVE approach rather resembles the TSVD image with a cutoff of 0.0001%, i.e. regularization is too small. Therefore, it is advisable to compare different settings of the regularization parameter and make the final choice based on the abovementioned considerations.
Cortical LORETA
Cortical LORETA is principally the same technique as LORETA, however, Cortical LORETA is not computed in a 3D volume, but on the cortical surface.
The cortical reconstruction in BESA Research fed from BESA MRI is a closed 2D surface with no boundaries and a very close approximation of the actual cortical form. It consists of an irregular triangulated grid.
The Laplace operator that is used for identifying a smooth solution in a threedimensional space is exchanged with a Laplace operator that runs on the twodimensional cortical surface.
There is a wide variety of 2D Laplace operators with different characteristics. The general form of the discrete Laplace operator is
[math]\Delta f\left( p_{i} \right) = \frac{1}{d_{i}}\sum_{j \in N(i)}^{}{w_{ij}\left\lbrack f\left( p_{i} \right)  f\left( p_{j} \right) \right\rbrack},[/math]
where p_{i} is the ith node of the triangular mesh, f(p_{i}) is the value of a function f defined on the cortical mesh at the node p_{i}, w_{ij} is the weight for the connection between the nodes p_{i} and p_{j} and d_{i} is a normalization factor for the ith row of the operator. Furthermore, N(i) is the set of indices corresponding to the direct (also called "1ring") neighbors of p_{i}.
BESA offers the choice of three Laplace operators with slightly different characteristics.
 Unweighted Graph Laplacian: This is the simplest operator. It takes into account only the adjacency of the nodes and not the geometry of the mesh:
[math]
w_{ij} =
\begin{cases}
1, & \text{if } p_{i} \text{ and } p_{j} \text{ are connected by an edge} \\
0, & \text{otherwise}
\end{cases}
[/math]
[math]d_{i} = 1[/math]
 Weighted Graph Laplacian: This operator is similar to the unweighted graph Laplacian but with different weights for the different connections. The connections between nearby nodes get larger weights than the connections between farther nodes:
[math]w_{ij} = \frac{1}{\operatorname{dist}\left( p_{i},p_{j} \right)}[/math]
[math]d_{i} = \sum_{j \in N(i)}^{} {\operatorname{dist}\left(p_{i}, p_{j} \right)}[/math]
where dist (p_{i} , p_{j}) is the distance between the nodes p_{i} and p_{j}.
 Geometric Laplacian with mixed area weights: This operator takes into account the angles in the corresponding triangles into account as well as the area around the nodes in order to determine the connection weights:
[math]w_{ij} = \frac{\cot\left( \alpha_{ij} \right) + \cot\left( \beta_{ij} \right)}{2}[/math]
[math]d_{i} = A_{\text{mixed}}[/math]
where α_{ij} and β_{ij} denote the two angles opposite to the edge (i , j) and A_{mixed} is either the Voronoi area, or 1/2 of the triangle area or 1/4 of the triangle area depending on the type of the triangle.
Regularization and other parameters:
 SVD cutoff: The regularization for the inverse operator as a percent of the largest singular value.
 Depth weighting: Turn depth weighting on or off.
 Laplacian type: Selection of Laplacian operators (see above).
Notes:
 Starting Cortical LORETA: Cortical LORETA can be started from the submenu Surface Image of the Image menu or from the Image Selection button.
 Please refer to Chapter “Regularization of distributed volume images” for important information on regularization of distributed inverses.
Cortical CLARA
Cortical CLARA is principally the same technique as CLARA, but Cortical CLARA is not computed in a 3D volume, but on the cortical surface. Instead of using a LORETA image as the basis for the iterative application, cortical CLARA uses cortical LORETA.
Regularization and other parameters:
 SVD cutoff: The regularization for the inverse operator as a percent of the largest singular value.
 Depth weighting: Turn depth weighting on or off.
 Laplacian type: Selection of Laplacian operators (see Cortical LORETA).
 No of iterations: Number of iterations for CLARA. The more iterations are used, the sparser becomes the solution.
 Automatic: The algorithm tries to determine the number of iterations automatically. The goodness of fit (GOF) is calculated after every iteration and if there is a big jump in the GOF then the algorithm will stop. If no jumps appear during the calculations then CLARA iterates until the specified number of iterations is reached.
 Regularize iterations: If one wants to use different regularization for the CLARA iterations than the value specified as "SVD cutoff", this option should be selected.
 Amount to clip from img (%): Cortical CLARA uses the solution from the previous iteration as an additional weighting matrix for the current iteration. That weighting matrix is constructed by cutting the "low" activity from the solution. This number specifies how much of the activity should be cut from the previous solution in order to construct the weighting matrix. This value is given as a percentage of the maximal activity. Default value is 10%.
Notes:
 Starting Cortical CLARA: Cortical CLARA can be started from the submenu Surface Image of the Image menu or from the Image Selection button.
 Please refer to Chapter “Regularization of distributed volume images” for important information on regularization of distributed inverses.
Cortex Inflation
The inflated cortex is a smoothened version of the individual cortical surface with minimal metric distortions (Fischl, B. et al. (1999). Cortical SurfaceBased Analysis: II: Inflation, Flattening, and a SurfaceBased Coordinate System. NeuroImage, 9(2), 195–207). Gyri and sulci are smoothened out. The original distances between each point on the cortex and its neighbors are, however, mostly preserved.
Cortical LORETA map overlaid on top of the inflated cortical surface.
A lighter gray color overlaid on top of the surface image indicates the location of a gyrus of the individual cortex surface, while a darker gray color indicates the location of a sulcus. The inflated cortical surface can be computed in BESA MRI 2.0. For more details please refer to the BESA MRI 2.0 help.
Surface Minimum Norm Image
Introduction
The minimum norm approach is a common method to estimate a distributed electrical current image in the brain at each time sample (Hämäläinen & Ilmoniemi 1984). The source activities of a large number of regional sources are computed. The sources are evenly distributed using 1500 standard locations 10% and 30% below the smoothed standard brain surface (when using the standard MRI) or using between 30004000 locations on the individual brain surface defined by the graywhitematter boundary.
Since the number of sources is much larger than the number of sensors in a minimum norm solution, the inverse problem is highly underdetermined and must be stabilized by a mathematical constraint, the minimum norm. Out of the many current distributions that can account for the recorded sensor data, the solution with the minimum L2 norm, i.e. the minimum total power of the current distribution is displayed in BESA Research.
First, the forward solution (leadfield matrix L) of all sources is calculated in the current head model. Then, the source activities S(t) of all source components are computed from the data matrix D(t) using an inverse regularized by the estimated noise covariance matrix:
[math]\mathrm{S}\left( t \right) = \mathrm{R} \cdot \mathrm{L}^{T} \cdot \left( \mathrm{L} \cdot \mathrm{R} \cdot \mathrm{L}^{T} + \mathrm{C}_N \right)^{1} \cdot \mathrm{D}(t)[/math]
Here, L is the leadfield matrix of the distributed regional source model, C_{N} denotes the noise correlation matrix in sensor space, and R is a weighting matrix in source space. R and C_{N} can be designed in different ways in order to optimize the minimum norm result. The total activity of each regional source is computed as the root mean square of the source activities S(t) of its 3 (MEG:2) components. This total source activity is transformed to a colorcoded image of the brain surface. (When the standard brain is used, two sources are assigned to each surface location, located 10% and 30% below the surface, respectively. The color that is displayed on the standard brain surface is the larger of the two corresponding source activities.)
Weighting options
The minimum norm current imaging techniques of BESA Research provide different weighting strategies. Two weighting approaches are available: Depth weighting and spatiotemporal approaches.
 Depth weighting: Without depth weighting, deep sources appear very smeared in a minimumnorm reconstruction. With depth weighting, both deep and superficial sources produce a similar, more focal result. If this weighting method is selected, the leadfield of each regional source is scaled with the largest singular value of the SVD (singular value decomposition) of the source's leadfield.
 Spatiotemporal weighting: Spatiotemporal weighting tries to assign large weight to sources that are assumed to be more likely to contribute to the recorded data.
 Subspace correlation after single source scan: This method divides the signal into a signal and a noise subspace. The correlation of the leadfield of a regional source i with the signal subspace (p_{i}) is computed to find out if the source location contributes to the measured data. The weighting matrix R becomes a diagonal matrix. Each of the three (MEG: 2) components of a regional source get the same weighting value p_{i}^{2}. This approach is based on the signal subspace correlation measure introduced by J.C. Mosher, R. M. Leahy (Recursive MUSIC: A Framework for EEG and MEG Source Localization, IEEE Trans. On Biomed. Eng. Vol. 45, No. 11, November 1998)
 Dale & Sereno 1993: In the approach of Dale and Sereno (J Cogn Neurosci, 1993, 5: 162176) a signal subspace needs not be defined. The correlation p_{i }of the leadfield of regional source i with the inverse of the data covariance matrix is computed along with the largest singular value λ_{max }of the data covariance matrix. The weighting matrix R is a diagonal matrix with weights: . Each of the three (MEG: 2) components of a regional source receives the same weighting value.
Noise regularization
Two methods to estimate the channel noise correlation matrix C_{N} are provided by the program:
 Use baseline: Select this option to estimate the noise from the userdefinable baseline. The signal is computed from the data at nonbaseline latencies.
 Use 15% lowest values: The baseline activity is computed from the data at those 15% of all displayed latencies that have the lowest global field power. The signal is computed from all displayed latencies.
In each case, the activity (noise or signal, respectively) is defined as rootmeansquare across all respective latencies for each channel.
The noise covariance matrix C_{N} is constructed as a diagonal matrix. The entries in the main diagonal are proportional to the noise activity of the individual channels (if selected) or are all equally proportional to the average noise activity over all channels. The noise covariance matrix C_{N} is then scaled such that the ratio of the Frobenius norms of the weighted leadfield projector matrix (LRL^{T}) and the noise covariance matrix C_{N} equals the SignaltoNoise ratio. This scaling can be multiplied by an additional factor (default=1) to sharpen (<1) or smoothen (>1) the minimum norm image.
Applying the Minimum Norm Image
The minimumnorm algorithm is started via the Surface minimum norm image dialog box, which is opened from the Image menu, or by typing the shortcut CtrlM: Please refer to Chapter “Surface Minimum Norm Tab” for more details.
As opposed to the other 3D images available from the Image menu, the surface minimum norm image is not computed on a volumetric grid, but rather for locations on the brain surface. Accordingly, the results of the minimum norm image are displayed superimposed to the brain surface mesh rather than to the volumetric MR image.
The figure below shows a minimum norm image computed from the file Examples\Epilepsy\Spikes\SpikesChild4_EEG+MEG_averaged.fsg. The EEG spike peak was imaged using the individual brain surface of the subject. A baseline from 300 to 70 ms was used. Minimum norm was computed with depth weighting, Spatiotemporal weighting according to Dale & Sereno 1993 and individual noise weighting with a noise scale factor of 0.01. The minimum norm image reveals the location of the spike generator in the close vicinity of the frontal lefthemispheric lesion in this subject.
Multiple Source Probe Scan (MSPS)
Introduction
The MSPS function provides a tool for the validation of a given solution. It is based on the following theoretical consideration: If the recorded EEG/MEG data has been modeled adequately, i.e. all active brain regions are represented by a source in the current solution, then any additional probe source added to the solution will not show any activity apart from noise. The only exception occurs if this probe source is placed in close vicinity to one of the sources in the current solution. In that case, the solution's source and the probe source will share the activity of the corresponding brain area. The MSPS applies these considerations by scanning the brain on a predefined grid with a regional probe added to the current solution. Grid extent and density can be specified in the Image settings. The power P of the probe source at location r in the signal interval is compared with the power of the probe source in a reference interval, defining a value q:
[math]\mathrm{q}\left( r \right) = \sqrt{\frac{\mathrm{P}\left( r \right)}{\mathrm{P}_{\text{ref}}\left( r \right)}}  1[/math]
MSPS can be computed on time domain or timefrequency domain data:
 In the time domain, q(r) is computed from the source waveform of the probe source. Here, P(r) is the mean power of the probe source at location r in the marked latency range, and P_{ref}(r) is the mean probe source power in the userdefinable baseline interval.
 In the timefrequency domain, an MSPS image can be computed from the complex cross spectral density matrices. By applying the inverse operator for a source configuration consisting of the current solution and the probe source, the power of the probe source can be computed for the target interval [P(r)] and the reference timefrequency interval [P_{ref}(r)]. In the resulting MSPS image, qvalues are shown in %, where q[%] = q*100.
The inverse operator used to determine the probe source power uses different regularization constants for the probe source and the sources in the current solution. The regularization constant of the sources in the current solution can be specified in the Image settings (default 4%). The regularization constant of the probe source is internally set to 0%.
Alternatively to the definition above, q can also be displayed in units of dB:
[math]\mathrm{q}\left\lbrack \text{dB} \right\rbrack = 10 \cdot \log_{10}\frac{\mathrm{P}\left( r \right)}{\mathrm{P}_{\text{ref}}\left( r \right)}[/math]
Values of q smaller than zero are not shown in the MSPS image.
According to the considerations above, an MSPS of a correct source model should optimally yield image maxima around the sources in the current solution only. If the MSPS image is blurred or shows maxima at locations different from the modeled sources, this indicates a nonsufficient or incorrect solution.
Applying the MSPS
This chapter illustrates the application of the Multiple Source Probe Scan. The figures are generated with data from file Examples/Epilepsy/Spikes/RolandicSpikeChild.fsg (300 : +200 ms, filtered from 3 Hz [forward] to 40 Hz [zerophase]).
Time domain versus timefrequency domain MSPS
The multiple source probe scan can be computed in the time domain or the timefrequency domain. The latter is possible only when timefrequency domain data is available for the current condition, i.e. if the condition has been created by starting a multiple source beamformer (MSBF) computation from the source coherence window. In this case, evoking the MSPS calculation from the Imaging button or from the Image menu will bring up the following dialog window that allows to choose between time or timefrequency MSPS. If only time domain data is available, this dialog window will not appear and MSPS will be computed in the time domain.
For a timefrequency domain MSPS, the target and the reference timefrequency interval have been specified already in the TimeFrequency window (see Chapter "How To Create Beamformer Images"). For a timedomain MSPS, the target and the reference epoch have to be specified in the Source Analysis window as described below.
Time domain MSPS
The timedomain MSPS image displays the ratio of the power of a regional probe source in the signal and the baseline interval. The currently set baseline is indicated by a horizontal line in the upper left corner of the channel box.
By default, BESA Research defines the prestimulus interval of the current data segment as baseline. The baseline should represent a latency range in which no eventrelated activity is present in the data. There are several possibilities to modify the baseline interval: by clicking on the horizontal line with the left mouse button or by using the corresponding entry in the Condition menu or Fit Interval popup menu.
Mark an interval to define the target epoch, i.e. the timeinterval for which the current solution is to be tested. Start the MSPS by selecting it from the Image selection button or from the Image menu to start the probe source scan. The Image menu can be evoked either from the menu bar or by rightclicking anywhere in the source analysis window. The 3D window opens and displays the scan result.
TimeResolved MSPS
If the MSPS has been computed on time domain data, the image can be shown separately for each latency in the selected interval. After the MSPS has been computed for the marked epoch, doubleclick anywhere within this epoch to display the ratio of the probe source magnitude at the selected latency and the mean probe source magnitude in the baseline. Scanning the latency range by moving the cursor (e.g. with the left and right arrow cursor keys) provides a timeresolved MSPS image.
Timeresolved MSPS images are not available if the MSPS has been computed on data in the timefrequency domain.
Notes:
 You can hide or redisplay the last computed image by selecting the corresponding entry in the Image menu.
 The current image can be exported to ASCII or BrainVoyager vmpformat from the Image menu.
 For scaling options, please refer to the scaling buttons popup menu .
 Parameters used for the MSPS calculations can be set in the General Settings tab of the Image Settings dialog box.
Source Sensitivity
Introduction
The 'Source sensitivity' function displays the sensitivity of the selected source in the current source model to activity in other brain regions. Sensitivity is defined as the fraction of power at the scanned brain location that is mapped onto the selected source.
To compute the source sensitivity, unit brain activity is modeled at different locations (probe source) throughout the brain. To this data, the current source model is applied to compute the source waveforms S_{CM} of all modeled sources:
[math]\mathrm{S}_{\text{CM}} = \mathrm{L}_{\text{CM}}^{1} \cdot \mathrm{L}_{\text{PS}}[/math]
Here L_{CM}^{1} is the regularized inverse operator for the current model, and L_{PS} is the leadfield of the regional probe source (dimension [Nx3] for EEG and [Nx2] for MEG, respectively, where N is the number of sensors). The source amplitude S_{SS }of the selected source in the model is a 3x3 (MEG: 2x2) submatrix of S_{CM} (if the selected source is a regional source) or a 1x3matrix (MEG: 1x2) (if the selected source is a dipole). The root mean square of the singular values of S_{CM} is defined as the source sensitivity.
The 3D source sensitivity image displays this value for all locations on a grid specified under Image/Settings. Grid density can be specified in the Image Settings.
Applying the Source Sensitivity Image
The Source Sensitivity image is evoked from the Image menu or by pressing the corresponding hot key (default: V).
This function is enabled only when a solution with an active selected source is present in the Source Analysis window. The source sensitivity image then displays the sensitivity of the selected source to activity in other brain regions.
Source Sensitivity image for the selected frontal source (green) in model 'High_Intensity_3RS.bsa' in folder 'Examples/ERP_Auditory_Intensity'. The data displayed is the '100dB' condition in file 'All_Subjects_cc.fsg'. The selected source is sensitive to activity in the frontal brain region (yellow/white), while it is not influenced by activity in the vicinity of the left and right auditory cortex areas, which are modeled by the red and blue source in the model (transparent/gray).
Notes:
 The sensitivity image is independent of the recorded sensor signals. It only depends on the current source model, the sensor configuration, the head model, and the regularization constant.
 If the regularization constant is set to zero, each source has a sensitivity of 100% to activity around its own location. With increasing regularization, the spatial filter becomes less focused, and the sensitivity of a source to activity at its location decreases.
 You can hide or redisplay the last computed image by selecting the corresponding entry in the Image menu.
 The current image can be exported to ASCII or BrainVoyager vmpformat from the Image menu.
SESAME
(requires Besa Research 7.0 or higher)
SESAME (Sequential SemiAnalytic MonteCarlo Estimation) is a Bayesian approach for estimating sources that uses MarkovChain MonteCarlo method for efficient computation of the probability distribution as described in Sommariva, S., & Sorrentino, A. "Sequential Monte Carlo samplers for semilinear inverse problems and application to magnetoencephalography." Inverse Problems 30.11 (2014): 114020.
It allows to automatically estimate simultaneously the number of dipoles, their locations and time courses requiring virtually no user input. The algorithm is divided in two blocks:
 The first block consists of a Monte Carlo sampling algorithm that produces, with an adaptive number of iterations, a set of samples representing the posterior distribution for the number of dipoles and the dipole locations.
 The second block estimates the source time courses, given the number of dipoles and the dipole locations.
The Monte Carlo algorithm in the first block works by letting a set of weighted samples evolve with each iteration. At each iteration, the samples (a multidipole state) approximates the nth element of a sequence of distributions p1, …, pN, that reaches the desired posterior distribution (pN = p(xy)). The sequence is built as pN = p(x) p(yx) α(n), such that α(1) = 0, α(N) = 1. The actual sequence of values of alpha is determined online. Dipole moments are estimated after the number of dipoles and the dipole locations have been estimated with the Monte Carlo procedure. This continues until a steady state is reached.
The SESAME image in BESA Research displays the final probability of source location along with an estimate for number of sources. Using the menu function Image / Export Image As... you have the option to save this SESAME image.
Notes:
 Grid spacing: Due to memory and computational limitations, it is recommended to use SESAME with a grid spacing of 5 mm or more.
 Fit Interval: SESAME requires a fit interval of more than 2 samples to start the computation.
 Computation time: Computation speed during SESAME calculation depends on the grid spacing (computation is faster with larger grid spacing) and number of channels.
Brain Atlas
(requires Besa Research 7.0 or higher)
Introduction
Brain atlas is a priori data that can be applied over any discrete or distributed source image displayed in the 3D window. It is a reference value that strongly depends on the selected brain atlas and should not be used as medical reference since individual brains may differ from the brain atlas.
Brain Atlases
In BESA Research the atlases listed below are provided. BESA is not the author of the atlases; please cite the appropriate publications if you use any of the atlases in your publication.
Brainnetome
This is one of the most modern brain probabilistic atlas where structural, functional, and connectivity information was used to perform cortical parcellation. It was introduce by Fan and colleagues (2016), and is still work in progress. The atlas was created using data from 40 healthy adults taking part in the Human Connectome Project. In March 2018, the atlas consists of 246 structures labeled independently for each hemisphere. In BESA we provide the max probability map with labeling. Please visit the Brainnetome webpage to see more details related to the indicated brain regions (i.e. behavioral domains, paradigm classes and regions connectivity).
AAL
Automated Anatomical Labeling atlas was created in 2002 by TzourioMazoyer and collegues (2002). It is the mostly used atlas nowadays. The atlas is based on the averaged brain of one subject (young male) who was scanned 27 times. The atlas resolution is 1 mm isometric. The brain sulci were drawn manually on every 2mm slice and then brain regions were automatically assigned. The atlas consists of 116 regions which are asymmetrical between hemispheres. The atlas is implemented as in the SPM12 toolbox.
Brodmann
The Brodmann map was created by Brodmann (1909). The brain regions were differentiated by cytoarchitecture of each cortical area using the Nissi method of cell staining. The digitalization of the original Brodmann map was performed by Damasio and Damasio (1989). The digitalized atlas consists of 44 fields that are symmetric between hemispheres. BESA used the atlas implementation as in Chris Roden’s MRICro software.
AAL2015
Automated Anatomical Labeling revision 2015. This is the updated AAL atlas. In comparison to the previous version (AAL) mainly the frontal lobe shows a higher degree of parcellation (Rolls, Joliot, and TzourioMazoyer 2015). The atlas is implemented as in the SPM12 toolbox.
Talairach
Atlas was created in 1988 by Talairach and Tournoux (1988) and it is based on the post mortem brain slices of a 60 year old right handed European female. It was created by drawing and matching regions with the Brodmann map. The atlas is available at 5 tissue levels, however we used only the volumetric gyrus level as it is the most known in neuroscience and is the most appropriate for EEG. The atlas consists of 55 regions that are symmetric between hemispheres. The native resolution of the atlas was 0.43x0.43x25 mm. Please note that the poor resolution in Z direction is a direct consequence atlas definition, and since it is a postmortem atlas it will not correctly match the brain template
(noticeable mainly on brain edges). The atlas digitalization was performed by Lancaster and colleagues (2000) resulting in a “golden standard” for neuroscience. The atlas was first implemented in a software called talairach daemon.
Visualization modes
Just Labels
Displayed are only Talairach Coordinates, the currently used brain atlas and the region name where the crosshair is placed. No atlas overlay will be visible on the 3D image.
brainCOLOR
All information is displayed as in “Just Labels” mode but also the atlas is visible as an overlay over the MRI. The coloring is performed using the algorithm introduced by Klein and colleagues (Klein et al. 2010). With this method of coloring the regions which are part of the same lobe are colored in a similar color but with different color shade. The shade is computed by the algorithm to make these regions visually differentiable from each other as much as possible.
Individual Color
In this mode the native brain atlas color is used if provided by the authors of the brain atlas (i.e. Yeo7). Where this was not available BESA autogenerated colors for the atlas using an approach similar to political map coloring. This approach aims to differentiate most regions that are adjacent to each other and no presumptions on lobes is applied.
Contour
Only region contours (borders between atlas regions) are drawn with blue color. This is the default mode in BESA Research.
References
 Brodmann, Korbinian. 1909. Vergleichende Lokalisationslehre Der Großhirnrinde. Leipzig: Barth. https://www.livivo.de/doc/437605.
 Damasio, Hanna, and Antonio R. Damasio. 1989. Lesion Analysis in Neuropsychology. Oxford University Press, USA.
 Fan, Lingzhong, Hai Li, Junjie Zhuo, Yu Zhang, Jiaojian Wang, Liangfu Chen, Zhengyi Yang, et al. 2016. “The Human Brainnetome Atlas: A New Brain Atlas Based on Connectional Architecture.” Cerebral Cortex 26 (8): 3508–26. https://doi.org/10.1093/cercor/bhw157.
 Klein, Arno, Andrew Worth, Jason Tourville, Bennett Landman, Tito Dal Canton, Satrajit S. Ghosh, and David Shattuck. 2010. “An Interactive Tool for Constructing Optimal Brain Colormaps.” https://mindboggle.info/braincolor/colormaps/index.html.
 Lancaster, Jack L., Marty G. Woldorff, Lawrence M. Parsons, Mario Liotti, Catarina S. Freitas, Lacy Rainey, Peter V. Kochunov, Dan Nickerson, Shawn A. Mikiten, and Peter T. Fox. 2000. “Automated Talairach Atlas Labels for Functional Brain Mapping.” Human Brain Mapping 10 (3): 120–131.
 Rolls, Edmund T., Marc Joliot, and Nathalie TzourioMazoyer. 2015. “Implementation of a New Parcellation of the Orbitofrontal Cortex in the Automated Anatomical Labeling Atlas.” NeuroImage 122 (November): 1–5. https://doi.org/10.1016/j.neuroimage.2015.07.075.
 Talairach, J, and P Tournoux. 1988. CoPlanar Stereotaxic Atlas of the Human Brain. 3Dimensional Proportional System: An Approach to Cerebral Imaging. Thieme.
 Thomas Yeo, B. T., F. M. Krienen, J. Sepulcre, M. R. Sabuncu, D. Lashkari, M. Hollinshead, J. L. Roffman, et al. 2011. “The Organization of the Human Cerebral Cortex Estimated by Intrinsic Functional Connectivity.” Journal of Neurophysiology 106 (3): 1125–65. https://doi.org/10.1152/jn.00338.2011.
 TzourioMazoyer, N., B. Landeau, D. Papathanassiou, F. Crivello, O. Etard, N. Delcroix, B. Mazoyer, and M. Joliot. 2002. “Automated Anatomical Labeling of Activations in SPM Using a Macroscopic Anatomical Parcellation of the MNI MRI SingleSubject Brain.” NeuroImage 15 (1): 273–89. https://doi.org/10.1006/nimg.2001.0978.
Review  

Source Analysis 

Integration with MRI and fMRI  
Source Coherence  
Export  
MATLAB Interface  
Special Topics 