Mapping electrical activity in the human neocortex - janeshdev/ggplot2 GitHub Wiki
The amount of light absorbed and scattered by neocortical (surface) brain tissue is altered by neuronal activity. Imaging of intrinsic optical signals (ImIOS), a technique for mapping these optical changes with an imaging detector, has the potential to be useful for both clinical and experimental investigations of the human neocortex. In particular, it may be possible in the future for neurosurgeons to use ImIOS to locate those regions of the neocortex that respond abnormally to electrical stimuli, and which therefore may be implicated in epileptic seizures. However, the usefulness of ImIOS for human studies is currently limited because the data are noisy and contaminated with respiration and heartbeat artifacts. In this case study we display results from a Bayesian, dynamic linear modeling (DLM) approach that appears to address these problems. The example shows that a smooth signal can be extracted from the noisy data and that the signal has the potential to help us find those regions of the brain that respond abnormally.
The image below is a human brain during surgery. Most of the image is gray matter overlain with a network of dark blood vessels. The brain is under a glass plate, for stability. The two light gray circles near the center are electrodes that are used to stimulate the brain. The dark gray circle just above them is an electrode for recording EEG signals, which are a direct measurement of the electrical activity of nearby neurons. Finally, several regions of the brain are numbered in red.
To learn which regions are hyperexcitable, and therefore possibly sites where epileptic seizures originate, the surgeon passes electrical current between the two light gray electrodes. ImIOS is based on what happens next. The electrical signal spreads across the neocortex through the firing of neurons. When neurons fire they require more oxygen, so they release molecules that signal nearby small arteries to dilate and bring more blood to the region. The blood carries oxygen to the neurons, but it also carries hemoglobin, which absorbs light. Digital photographs can record which areas become darker and are, therefore, more neuronally active.
The data we work with come from a series of digital photographs acquired several times per second, starting about 10 seconds prior to stimulation and lasting about two minutes. For each region, we get a time series y_1, ..., y_T where y_t is the average percentage change from baseline of all the pixels in a region. T is typically between about 500 and 1000.
The next figure shows the first 50 seconds of data collected from one region. The rapid oscillations, about every 2-3 time steps, are due to heartbeat. The peaks at around t = 3, 8, 14, 19, ..., 40, 46 seconds are due to respiration.
To filter out the effects of heartbeat and respiration and extract the smooth, underlying signal, we use a Normal, linear state-space model having a locally linear component for the signal and quasi-cyclical components for heartbeat and respiration. The figures below are typical of the results. The top figure shows the raw data and the estimated signal for two regions. The bottom figure shows the estimated heartbeat (red) and respiration (green) components and an additional (blue), quasi-cyclic component for which we do not know the physiological explanation.
The surgeon typically stimulates the brain several times, with currents of different strength (amps). The next figure shows data from four regions, when stimulated once at high intensity "a" and once at low intensity "c". The stimulus was applied for about four seconds beginning around t=20, indicated by the gray rectangle in each facet. The dark trace is from the EEG electrode, which records electrical activity of the neurons near the stimulus. For replicate "c", the EEG is very active when the stimulus is on but returns to normal almost immediately. For replicate "a", the EEG remains active for many seconds after the stimulus and dies off slowly, so-called "after-discharge activity". For surgery, it is important to learn which regions exhibit after-discharge activity and how strong the current must be to induce it.
ImIOS can help if the underlying signal has different characteristics according to whether there is after discharge activity. Some differences can be seen in the previous plot of y_t vs. time: the curves from replicate "a" appear to decrease more gradually and have a more rounded minimum than those from "c". Another way to view the signal is in a phase portrait --- a plot of slope versus level, as shown next. In this figure, time is coded by color, to match the previous figure. The phase portraits all start in red near the origin, transition from yellow to green during stimulus, and finally return to a resting state in blue and magenta. Here, some the curves from "a" take a turn to the left during stimulus. It remains for further investigation to determine whether these aspects of ImIOS can reliably indicate the presence of after-discharge activity.
The code for the first several plots is straightforward, so we show the code just for the last two.
Data are in three data frames.
p4.535.m: ImIOS data. contains time, rep, region, value
fitted.simple.p4.535.m: model output. contains time, rep, region, param, value
eeg.m: EEG data. contains time, rep, eeg
The first plot consists of four layers. The first is geom_polygon, which creates the gray
rectangle. It has to come first so the others will plot over it. Second is geom_point, which
plots the raw ImIOS data. Third is geom_line which adds the smooth curve estimated by the
model. Fourth is another geom_line to add the eeg trace. The eeg trace begins 18.78 seconds after
the ImIOS data, so 18.78 is added to its x-coordinate. We displaced it upward by adding 7 to its
y-coordinate so it does not obscure the data. The plot begins with geom="blank", to set up
some default values for subsequent layers. We mapped color to time, with
scale_color_gradientn, so we can match times between this plot and the next.
p1 <- ( qplot ( time, value, geom = "blank",
data = subset ( p4.535.m, region %in% 1:4 & time <= 100 ) )
+ geom_polygon ( data = data.frame ( x = c ( 20.72, 20.72, 24.67, 24.67 ),
y = c ( -15, 15, 15, -15 )
),
mapping = aes ( x = x, y = y ),
fill = I("gray86")
)
+ geom_point ( data = subset ( p4.535.m, region %in% 1:4 & time <= 100 ),
mapping = aes ( x = time, y = value, color = time ),
shape = 1, alpha = 1/5)
+ geom_line ( data = subset ( fitted.simple.p4.535.m,
param == "level" & time <= 100 ),
mapping = aes ( x = time, y = value, color = time ) )
+ geom_line ( data = eeg.m, mapping = aes ( time+18.78, eeg+7 ) )
)
( p2 <- p1 + scale_y_continuous ( expression ( y[t] ) )
+ scale_x_continuous ( "time (sec)" )
+ facet_grid ( region ~ rep )
+ theme_bw()
+ scale_color_gradientn ( colours = rainbow(12) )
)
For the second plot, first we recast fitted.simple.p4.535.m
in a form more convenient for the plot. And second, we mapped color to time, with a
scale identical to the first plot, so we can see which parts of the phase portraits
correspond to which parts of the time series.
foo <- cast ( fitted.simple.p4.535.m, time + rep + region ~ param )
jpeg ( "rainbow2.jpg", quality=100 )
( qplot ( level, slope, data = subset ( foo, time <= 100 ),
geom="path", color = time )
+ facet_grid ( region ~ rep )
+ theme_bw()
+ scale_color_gradientn ( colours = rainbow(12) )
)web page created by
Michael Lavine
Department of Mathematics and Statistics
UMass Amherst
Amherst, MA, USA
data acquired by
Michael Haglund and Daryl Hochman
Department of Neurosurgery
Duke University
Durham, NC, USA
Data were acquired under an IRB protocol which, to protect the privacy and rights of human subjects, does not permit posting on a web site.