The purpose of this page is just to serve as todo or scratch pad for the development project and to list and share some ideas.

After making changes to the code and/or documentation, this page should remain on the wiki as a reminder of what was done and how it was done. However, there is no guarantee that this page is updated in the end to reflect the final state of the project

So chances are that this page is considerably outdated and irrelevant. The notes here might not reflect the current state of the code, and you should not use this as serious documentation.

Background on the dataset

The purpose of the study in which the data that is demonstrated here is to use EEG with Optogenetic Stimulation (opto-EEG) for studying sensorymotor integration.

Sensorymotor integration is a neurological process dealing with somatosensation inputs into motor outputs in an effective way. Numerous studies have shown the direct signal pathways within somatosensor-motor circuits, however little is known how the tactile sensations with different stimulation frequency are processed from somatosensory to motor cortex in order to react appropriately to different types of stimulation. In present study, we mimicked the peripheral sensation by direct stimulation of S1, S2, M1 and sensory thalamus using optogenetic method and concurrently recorded the frequency dependent responses (1, 10, 20, 30, 40 and 50 Hz) with a depth electrode in the region of the optode (i.e. S1, S2, M1 and thalamus), and recording the EEG using the surface micro electrode array.

Subject and Experiment

In the study applied here, we used the transgenic mice (Thy1-ChR2-EYFP, B6 mice, 12 weeks; body weight 27 g, malehttp://jaxmice.jax.org/strain/007615.html). For implantation of electrodes on the skull and sensory thalamus (VPM; 1.82 mm posterior, 1.5 mm lateral and 3.7 mm ventral to bregma), mice were anesthetized (ketamine/xylazine cocktail, 120 and 6 mg/kg, respectively) and then positioned in a stereotaxic apparatus. Location of electrode placement is going to describe in section ?. We allocated two electrodes in the most anterior region as reference and ground electrodes. The brain signals are recorded by both high-density electrode array (EEG-38 channels + ground and reference-2 channels) and local field potential (single channel). The EEGs were acquired with an analog amplifier (Synamp, Neuroscan, USA) with a sampling frequency of 2000 Hz.

Preprocessing

define trials

Using the FieldTrip function ft_definetrial you can define the pieces of data that will be read in for preprocessing. Trials are defined by their begin and end sample in the data file and each trial has an offset that defines where the relative t=0 point (usually the point of the optogenetic stimulus-trigger) is for that trial.
The ft_definetrial and ft_preprocessing functions require the original EEG dataset acquiring from Synamp system, which is available at this link <cnt file link>.

Actually, example dataset has not digital trigger information. To mark stimulation timing in the file, they use analog input (41th channel) as trigger information. Since FieldTrip can offer to support customized function by using cfg.trialfun, this tutorial shows how to make trial-based dataset.
This results in a cfg.trl = 'mousetrialfun' in which the beginning, the trigger offset and the end of each trial relative to the beginning of the raw data is defined.

You will see that the trial function produces a figure, which can be used to check whether the threshold was at the appropriate level and whether the onset of stimulation is properly detected.

insert figure (1 data loading)

The segments of interest are specified by their begin- and end-sample and by the offset that specifies the timing relative to the data segment. The offset is zero here indicating that the first sample is interpreted as time t=0.

preprocessing

After the segments of interest have been added to the configuration structure, we read them from the raw data file into memory. At this spet we also add the other preprocessing options to the configuration, such as the filter settings.

In the data.sampleinfo you can find the begin and endsample of each trial. The offset has been used to make an individual time axis for each trial.

We can simply plot all channels in a single trial using standard MATLAB code:

plot(data.time{1}, data.trial{1});
legend(data.label);
grid on

insert figure (2 single trial plot )

checking for artifacts

Using ft_databrowser we can do a quick visual inspection of the data and check whether there are artefacts.

cfg = [];
cfg = ft_databrowser(cfg, data);

insert figure (3 datatrial browser)

The HL1 channel contains the analog representation of the stimulus, which is of much larger amplitude than all other channels. We can exclude it in the guy (click “channel” button) or using the following:

We see that some channels show a much lower noise level than others (e.g. AF8). This might be due to differences in electrode-skull impedance.

insert figure (4 check the impedance and report it here)

rereferencing

The two most frontal electrodes in the grid were used as ground and reference during acquisition. These two electrodes are placed just anterior of the Fp1 and the Fp2 electrode.

This section shows a technique known as common average referencing (CAR) to generate a simple re-referencing method for array-type electrode recordings. CAR is a computationally simple technique. CAR is commonly used in EEG, where it is necessary to identify small signal sources in very noisy recordings.

We hypothesize that CAR will improve neural recording quality with respect to 38 ch electrode references.

Channel level analysis - ERPs

making a channel layout and deal with animal size

In order to make layout of high-density electrode array (HD-array) on mouse skull, FieldTrip support layout files ('mouse_layout.mat') that gives you exact control of the 2-D position of the sensors for topoplotting, and of the per-channel local coordinate axes for the multiplotting. The mat-file containing various variables to draw outlines of the headshape within the name “lay”.
Let us give you the information how to make customized layout for HD-array by using FieldTrip.
This figure indicates an electrode arrangement of HD-array (this example photo is taken from Lee et al. (2011) Journal of Visualized Experiments (video, pdf).

You can specify cfg.image in ft_prepare_layout and subsequently click on the location of each electrode. After specifying each electrode location, you'll be asked to specify the outlines of the head (i.e. the circle around the head, the eyes, the nose and ears (dash lines) and optionally some lines representing other important landmarks: bold lines, bregma, and lambda) and to specify the mask for the topographic interpolation (green dash line).
Find the bregma and lambda points on the skull, and record them with the stereotaxic ruler. The bregma and lambda points are located in the same anterio-posterior axis.
The anterio-posterior axis of high-density electrode array (HD-array) matches the line between bregma and lambda points. The bregma point (0, 0) should be located at the middle of the 4th layer of the anterior of HD-array because we follow the Paxinos coordinate system. We implemented a mouse layout on fixed length (4.2 mm) of between bregma and lambda.

After creating the layout, you should manually assign the correct name of the channel labels in the lay.label cell-array. Channel arrangement in lay.label is followed click sequence during ft_prepare_layout process. You can use ft_layoutplot for a visual inspection of the complete layout.

cfg = [];
cfg.layout = lay;
ft_layoutplot(cfg);

insert figure (6 layout plot)

If you think that some outlines (nose, eyes, head, whiskers and ears) are not necessary, you can pass without assigning any outline. Next step is zero point calibration for the bregma point. If you, however, use customized layout for single subject, you don’t need to carry out next step.

In the aspect of mouse anatomy, the bregma point is located in the middle of the 4th layer of the anterior of HD-array that it should be set (0, 0) because we follow the Paxinos coordinate system. To calibrate layout position, we just subtract the value of bregma on acquired layout from previous step.

HD-array is not scalable material. Each mouse, however, has different head size depending on strain, age, weight and gender. To circumvent this problem, we set to distance of bregma and lambda as the reference scale (4.2 mm). This distance could determine the level of head size with respect to individual subject. To get scalable mouse layout easily, we made a simple function ( ft_mouse_headsize_calibration ). If you use “ft_mouse_headsize_calibration” function, you can get calibrated head layout depending on the distance of bregma and lambda. The value of this function, cfg.B2Ldistance, is only affect size of outlines. This keeps the electrodes and the mask the same.
Take an example script as below:

The current standards to deal with differences in mouse brain size are very comparable to those adopted in the Talairach-Tournoux anatomical atlas of the human brain. For human EEG it is easy to make use of a template and to compare or average subjects, since EEG recording caps scale along with the size of the human head. E.g. at the Donders we have caps in different sizes, ranging from 52 to 60 cm head circumference in 2cm steps, thereby accommodating approximately a 15% variance in head circumference.

For the mouse the electrode grids are however of fixed dimensions, but the head sizes still differ.

For example, let's think about a case of targetting CA1 hippocampus which is (AP, ML, DV) = -2, 1.5, -2 with respect to the bregma according to the mouse atlas. Each idividual mouse has different brain size. What we do is to measure the length between bregma to lambda, and if the length is 4.2 mm, we multiply the target distance by 4.2/3.9 ~= 1.077. Hence the steretaxic target for CA1 would be (AP, ML, DV) = -2.15, 1.62, -2.15.
Since the tolerance range of the stereotaxic is 0.1 mm, the actual target would be -2.1, 1.6, -2.1.
I observed that some group rescale only anterior posterior (AP). Hence according to them the target coordinate would be -2.1, 1.5, -2. Unlike human, mouse has a secondary confirmation procedure, which is called histology. Hence, we can confirm whether we hit the target or not via histological procedure after recordings.
For your information, many neuroscientists agree that this rescaling method can work as a rule of thumb but individual difference between animals cannot be neglected, emphasizing the post-hoc histological confirmation.

The anatomical structure of mouse brain due not vary dramatically if the weight and/or the length between bregma and lambda are in the similar range, in a sense of a few hundreds micrometer. For example, any structure with a dimension ranging larger than a couple of hundreds um, we seldom fail to hit this structure. However, when we aim for a structure with a sharp shape with dimension smaller than a couple of hundreds um, we half fail to hit that structure. But if this structure is on the cortex or near, we seldome fail.

compute and visualize the ERPs

The ERP that is time locked to stimulus onset is computed by averaging the data segments over all trials.

Using the figure legend, we see that channel VPM shows a large response to the stimulation. This makes sense, because channel VPM corresponds to the depth electrode that is inserted along with the optode, i.e. it records from the stimulated site.

Rather than plotting all ERPs on top of each other, we can also plot them according to the channel layout.

You can use the zoom in and zoom out buttons of MATLAB. Furthermore, you can make a selection of channels and click in the selection, which causes them to be averaged and displayed in a single plot. In the single plot, you can make a selection of time (i.e. latency), which is subsequently averaged as an interpolated topographic distribution of the potential.

Channel level analysis - TFRs

In this section, we will describe how to calculate time-frequency representations using Hanning window.

visualizing TFRs

To visualize the power changes, a normalization with respect to a baseline interval will be performed.
This method, “relchange”, means (active period-baseline)/baseline. Note that the “relchange” is expressed as a ratio of subtracting, for each frequency.
If you want to visualize TFR on the head outline of mouse, “cfg.layout=mouse_layout.mat” have to be declared.

This gives a figure that shows the origin (white sphere), the x-, y-, and z-axes. We notice that the origin is not at Bregma, nor at the interaural point. This is inconsistent with the desired coordinate system, hence we have to realign the anatomical MRI.

Coregistration using high-level graphical function

We can use the ft_volumerealign function to coregister the anatomical MRI to the desired coordinate system. There are multiple coordinate systems in which the anatomy of the brain and head can be described, but here we want to use the Paxinos-Franklin cordite system which uses Bregma and Lambda as anatomical landmarks (or fiducials).

The anatomical MRI is displayed in three orthogonal plots. You have to visually identify the fiducials and press “b” for bregma, “l” for lambda and “z” for a midsagittal point. Using “f” you can toggle the fiducials on and off. Once you are happy with their placement, you press “q” and the realigned mri is returned

It is useful to also explicitly specify the coordinate system in the anatomical MRI. It is used in ft_sourceplot and various other functions to check whether various geometrical objects are expressed in the same coordinate system.

mri_realigned2.coordsys = 'paxinos';

Correcting the units of the anatomical MRI

Finally, we notice is that the units are off by a factor 10x. In the figure of ft_determine_coordsys, the sphere at the origin, each thick line segment is by default 1 cm, and the axes are 15 cm long, as that is appropriate for judging the human anatomy. We can also plot it with

The units can be fixed by

mri_realigned.unit = 'mm';

After coregistering the MRI with the Paxinos coordinate system, it is convenient to reslice it, i.e. to interpolate the greyscale values on a 3-D grid that is nicely aligned with the cardinal axes.

Since the original anatomical and labeled MRI are expressed in the same coordinate system, we can apply the same transformation to the atlas. Again, we also fix the units and specify the coordinate system.

You see that the anatomical labels are not what you would expect. This is due to the anatomically labeled MRI file not specifying the correct labels. We can fix this and manually assign the labels from this paper.

The skull shows up with integer value 1, the brain (which is fully enclosed in the skull) shows up with integer value 2.

insert figure (20 image segmentation)

make mesh of brain-skull and skull-skin boundary

This part describes how to make the volume conduction model for the source localization.
In Previous section we did volume segmentation to extract brain and skull. Since we use In vitro MRI without skull, virtual skull was made by image dilation method from brain surface.
To set up mesh objects, we use the ft_prepare_mesh to get the triangulated meshes for skull and brain.

We've been set limited number of vertices in ft_prepare_mesh.
Because OpenMEEG could not support to bigger volume above around 1500 vertices.

make BEM volume conduction model

After making up volume objects, we perform the ft_prepare_headmodel for assigning the electrical property of volumes. From the literature in human study, the brain conductivity ranges from 0.12-0.48/Ω.m [1-3], and the human skull from 0.006-0.015/Ω.m [4,5] or even higher as 0.032-0.080/Ω.m [5]. According to another studies for the conductivity ratio between skull and brain, they reported numerical value with large variation the ranges from 25 to 80 times [6]. It is hard to specify the brain-to-skull conductivity ratio from these values with such large variations. In this step, we just assign conductivities applying to 80 times ratio between skull and brain.

OpenMEEG is an external package that solves the forward problems. It implements a Boundary Element Method (BEM) and provides accurate solutions when dealing with realistic head models. As we mentioned above segmentation section, we are computing the two volume layers (skull and brain).

Specification of electrodes

Now, we have the volume conduction model (vol) from MRI and the channel information (loc) for array electrode within different coordinate.
This channel information is real distances of each electrode from middle of 4th layer (bregma point) in high-density electrode array.

Even after processing the co-registration, there could be have some gaps between electrode and skull surface. However, in the case of EEG in FieldTrip toolbox, the 3D electrode position is projected orthogonally onto the skull surface automatically.

make leadfield matrix

The final procedure of the forward problem is to generate a leadfield that representing the linear relation between sourcespace and measurements (Gain matrix).
By using ft_prepare_leadfield we can get the matrices with respect to aligned electrode position and BEM meshes.

Calculating the cross spectral density matrix

The beamforming technique is based on the spatial filter. The DICS spatial filter is derived from the frequency counterpart of the covariance matrix (the cross-spectral density matrix). This matrix contains the cross-spectral densities for all electrode combinations and is computed from the Fourier transformed data of the single trials. It is given as output when cfg.output = 'powandcsd'. The frequency of interest is 10 Hz and the smoothing window is +/-4 Hz:

Before applying cross-spectral density matrix, the redefine of time interval and the channel selection will be carry out like below.

Source reconstruction

Using the covariance matrices and the leadfield matrices a spatial filtering is calculated and estimated the dipole intensity for each grid point. By applying the filter to the Fourier transformed data we can then estimate the power for neural activity by optogenetic stimuli (10 Hz). This results in a power estimate for each grid point. To get normalized index for the neural activity we have to do the spatial filtering both pre-stimulus and post-stimulus.

To calculate neural activity index is the same like below equation. The function ft_sourceinterpolate interpolates the source reconstructed activity or a statistical distribution onto the voxels or vertices of an anatomical description of the brain (MRI with atlas).

Dealing with differences in size

Compared to human EEG, in the mouse EEG recordings the electrode grid does not scale along with the size of the animal. Consequently, the position of the electrodes relative to the brain depends on the head size. This affects the topographic plotting of channel level data, the comparison (group stats) of channel-level data from multiple animals, and the electrode positions and volume conduction model used for source reconstruction.

The principled solution to this is that the experimentally measured lambda-bregma distance (or some other measure of head size) is used to scale the background image in the channel layout and to scale the volume conduction model of the head.

demonstrate how to scale the outline and mask in the layout

An more pragmatic solution is to keep the head size the same, but rather to inverse-scale the electrodes and keep the head size constant. Effectively the result is the same, but it is easier to manage. Furthermore, it makes the source-level results directly comparable over animals.

demonstrate how to scale electrode positions, box width and height in the layout

demonstrate how to scale electrode positions for the source reconstruction