Walkthrough

This page is a tutorial on how to use BrainTrak on an existing data set. For information on how to prepare raw EEG data, see here.

Fitting a single spectrum

Our first task is to fit a single EEG spectrum. For this purpose, a set of example EEG data is provided in the file demo.mat. To fit and plot this example, use the following commands:

load demo.mat
a = bt.fit(bt.model.full,f,s_single)
a.plot()

This uses the function bt.fit(), which runs BrainTrak on a single EEG spectrum. The model object bt.model.full selects which model to use for fitting. The remaining arguments correspond to the EEG spectrum itself (which you can plot using loglog(f,s_single)). The output of the fit is a bt.feather object, which has a method plot() that creates a bt.viewer whose role is to manage a plot window like the one below:

The reason for bt.viewer being a class rather than just a plotting function is that the viewer contains a lot of additional functionality, which will be demonstrated a little further on in this walkthrough.

To fit a different model, you simply specify a different model object to the fit command. For example, to fit the reduced model, you can use:

b = bt.fit(bt.model.reduced,f,s_single)
b.plot()

Notice how running b.plot() creates a new, slightly different figure window. On the left, the spectrum and the XYZ plot are essentially the same as before. However, the marginal distributions on the right now correspond to the parameters from the reduced model. These plots are produced automatically based on the specification of the model parameters in the model object.

A common task is to compare the fit from two different models. Rather than using bt.feather.plot(), we can create a bt.viewer object directly. This allows us to create the object with two feather objects, rather than one. To superimpose the plots, use:

bt.viewer(a,b)

Now, the original viewer is displayed, but the spectrum plot has an additional spectrum corresponding to the fit in the second feather object.

Sometimes you might want to change the parameter values manually to see what effect they have on the fit. BrainTrak provides a 'widget' very much like the one in the corticothalamic-model repository, that allows you do interactively change the parameters. You can use it by calling the interactive_fit() method of the bt.feather object:

a.interactive_fit

This will open a widget like the one below:

You can do the same for the reduced model:

b.interactive_fit

Notice how the two widgets have different sliders corresponding to the particular parameters used by each of the models.

You can run the routine for longer by specifying the number of points. This can be an exact number of points, or a string corresponding to time in seconds. On my computer, 20 seconds yields 39000 points in the chain

c(1) = bt.fit(bt.model.reduced,f,s_single,'20s')

Or you can specify the number of points

c(2) = bt.fit(bt.model.reduced,f,s_single,1e4)

My computer can generate 10000 steps in about 5 seconds. Notice how feather objects can be stored in an array (no need for a cell array). In general, for publication and offline analysis it is more valid to specify a fixed number of points for all fits. However, for general investigation and real time fitting, it can be more useful to fix the total time that the routine will run for.

If you open a parallel pool, BrainTrak will automatically run in parallel. This will generally mean that you can simulate more points in the same amount of time. If you specify a fixed simulation time, you should get more points in your output. If you simulate a fixed number of points, the fit should be faster.

parpool
c(3) = bt.fit(bt.model.reduced,f,s_single,'20s')
c(4) = bt.fit(bt.model.reduced,f,s_single,1e4)

On my computer, the 20s run now gives me 140000 points, and the 10000 point run takes 3 seconds. You can use the array of feather objects to examine these quantities across the runs

arrayfun(@(x) x.fit_data.fit_time,c) % Print out all of the simulation run times
arrayfun(@(x) x.fit_data.chain_length,c) % Print out all of the chain lengths

Fitting a trajectory

The next major use for BrainTrak is in tracking parameters over time. This functionality is provided by bt.fit_track(). You can fit a trajectory for the demo data as shown below:

f = bt.fit_track(bt.model.full,'demo',1,1:10,'10s')
f.plot()

The arguments to bt.fit_track specify:

• The model (bt.model.full)
• The name of the data set ('demo')
• The subject index (1)
• The time indexes to fit (1:10 meaning the first 10 spectra in the data)
• The chain length (here, 10 seconds of computation)

bt.fit_track loads the data using bt.core.load_subject. Therefore, this file needs to know how to read your data. In particular, it needs to

• Know where the .mat files are stored
• Know which electrodes to use for your data
• Perform any preprocessing on the data e.g. truncating parts of the recording

Having run bt.fit_track, the output is a bt.feather object, the same as before. Thus you can generate a plot using f.plot().

This plot is very much like the ones from before, but now there is an additional plot along the bottom. This plot shows the goodness-of-fit over time. If you left click on the plot, it will start an animation that steps through each of the fits. The position along the animation is shown as a black vertical line. Left click again to stop playback. If you right click on the plot, the vertical bar will follow the mouse pointer. You can move your mouse left or right to scroll through the fit. Right click again to stop scrolling.

Given a feather object, you can create a new one with a subset of the data using bt.feather.subrange(). For example, you can retain and plot only the first 5 fits:

f2 = f.subrange(1:5)
f2.plot()

One of the uses of bt.viewer is that bt.feather retains a reference to the viewer that it created. The bt.feather.plot function also accepts an argument specifying the index of the fit to plot. If you run

f.plot(1)

it will open a viewer that shows only the first fit (with no navigation bar at the bottom). If you then run

f.plot(2)

then the fit will change in place to show the second fit. That is, the data will be updated without opening a new figure window.

You can plot the time courses of each of the fitted parameters using

f.timecourse

Again, the parameters that are plotted are read from the bt.model object, so if you use a model with different parameters, the timecourse figure will automatically show the parameters used by your selected model.

Another plot that can be useful to examine is a scatter plot in XYZ. You can make this using

f.point_cloud

The points above are color-coded according to the sleep stage, if sleep stage information is present in the data set - see here for more information. Blobs can be rendered by computing the 3D convex hull of the points.

f.blobs

These blobs are automatically computed independently for each of the sleep stages. That is, the points are grouped by sleep stage/color when computing the convex hulls. You can also plot a 3D 'histogram' by binning the data in XYZ, where the density of the cloud corresponds to the number of points in the bin

f.clouds

For all of these plots, the color scheme is provided in +bt_utils/state_cdata.m.

Spatial variations

For this next stage, add the folder 'data_examples' to your path. This folder contains wrapper scripts to fit multielectrode data from Brain Resource. an example of this data is contained in demo_spatial.mat. Note that unlike bt.fit_track, the function fit_br_spatial calls load_br_data which does equivalent preprocessing to bt.core.load_subject except customized to the BR multielectrode data. We can fit this data using

f = fit_br_spatial(bt.model.spatial_express_t0,'demo','30s')

Note that this takes quite a bit longer to run than the single electrode fits. Again, the output is a bt.feather object, so we can use the plot() method to display the plot.

f.plot()

Because the fit contains multielectrode data, all of the electrodes are automatically displayed. The fit can also be displayed using the physical electrode positions. To do this, use the head_plot() method of the feather object:

f.head_plot()

You can click on the electrodes to change which spectrum and fit are displayed. The spatial distribution of EEG power is shown on the left, at a frequency indicate by the vertical bar in the spectrum plot. You can click on the spectrum to change the frequency whose power distribution is plotted. The power distribution on the left corresponds to blue spectrum on the right - thus it shows the experimental power distribution. If you use

f.head_plot(true)

you can see that the blue and red traces are swapped. Now, the power distribution on the left side corresponds to the model predictions.

Next steps

• Want to know how the initial conditions in the standard corticothalamic model are obtained? Learn more about database fitting here
• Want to make a new model? First read about how models are implemented, then walk through creating a new one here
• Want to use BrainTrak with a new data set? First read about how data is processed for BrainTrak, then walk through an example here
• Want to generate new analyses and plots? Learn about feather objects and how BrainTrak outputs are stored and manipulated here
• Want to learn more about the fitting algorithm? Detailed information about how the adaptive MCMC system is implemented can be found here