Introduction to EZFragility package

Key definition of the neural fragility marker

Drug-resistant epilepsy (DRE) is a significant clinical challenge. Effective treatment requires the synthesis of large amounts of personalized clinical data with the goal of identifying the region within a patient’s brain from which seizures originate (the epileptogenic zone, or EZ). Current standard of care involves multidisciplinary group meetings and clinician consensus after reviewing patient data. This lack of a standardized methodology has led to variable surgical outcomes, measured by seizure freedom/reduction after surgery. Computational methods to identify the EZ have been developed to improve outcomes. The neural fragility implemented is the R package is a promising biomarker of the epileptogenic zone. We provide this R package so that the method can be tested on extensive data for research.

The hypothesis of the neural fragility is that small changes in connection strengths at seizure onset zone nodes (SOZ) cause an imbalance in inhibitory and excitatory connectivity between brain regions. Then if the SOZ is perturbed, over excitation can occur manifesting in a seizure. The neural fragility marker propose a metric to evaluate when connection strength changes cause imbalance and on which electrodes.To create the fragility marker, a personalized dynamic model of the brain network from observed iEEG signals is computed in a discrete number of time windows. Then the neural fragility measures the degree network nodes are imbalanced by systematically computing the minimum perturbation required for each electrode’s connections to produce instability for the entire network.

Load small ictal ieeg data for patient PT01

This example data corresponds to the first seizure of patient PT01 from the Fragility Data Set. pt01EcoG is a matrix with 3000 rows (time points) and 84 columns (electrodes). The acquisition frequency is 1000Hz and the time around seizure onset is [-1:2]s. The rows names are the electrodes names. The column names is the time relative to seizure onset in (s). For more information about the data, you can use the ?pt01EcoG command.

The patient underwent a surgical treatment and was seizure-free after the surgery. The electrodes located in the SOZ and the surgically removed area are marked in sozIndex attribute.

data("pt01EcoG")

## create an epoch object
epoch <- Epoch(pt01EcoG)
epoch
#> 84 x 3000 matrix of doubles: 
#> 
#>             [,1]        [,2]        [,3] ...    [,3000]
#> G1     16650.233   25734.890   37898.546 ...  96182.263
#> G2   -225911.037 -226337.719 -226012.026 ... 108857.517
#> G3   -230090.583 -235515.210 -245789.565 ... 149984.457
#> ...          ...         ...         ... ...        ...
#> SLT4  109452.929  103181.323  101214.372 ...  92960.977
#> 
#> Time range: -1 to 1.999
#> Use $ to access its methods. see "?`Epoch-method`"

The function visuIEEGData can visualize the iEEG data.

visuIEEGData(epoch  = epoch)

The Epoch object can be subset using the [ operator.

# The electrode names corresponding to the site of the patient's surgery
sozNames <- attr(pt01EcoG, "sozNames")

## Show the electrodes that are marked as SOZ and additional 4 electrodes
display <- c(sozNames, "MLT1", "MLT2", "MLT3", "MLT4")
visuIEEGData(epoch  = epoch[display])

## Equivalent to: 
## visuIEEGData(epoch  = epoch[display, ])

## constrain to the first 100 time points
visuIEEGData(epoch  = epoch[display, 1:100])

You can use truncateTime to specify the time range

epochClipped <- truncateTime(epoch, from = -1, to = 0)

visuIEEGData(epoch  = epochClipped)

Compute Fragility Matrix

The function calcAdjFrag computes the fragility matrix for the given iEEG time series. The simplest form is calcAdjFrag(ieegts, window, step). However, computing Fragility Matrix can take a bit while. To speed up the computation, we provided a parallel example below. The parallel feature is built upon the foreach package. You need to have a parallel backend registered in your R session beforehand. For example, you can use the doSNOW package to register a parallel backend. The result pt01Frag (data("pt01Frag")) from this code snippet corresponds to the data saved in the package. The function calcAdjFrag allows the user to tune the parameters window, step and nSearch (see help ?calcAdjFrag).

## Register a SNOW parallel backend with 4 workers
library(doSNOW)
#> Warning: package 'doSNOW' was built under R version 4.2.3
#> Loading required package: foreach
#> Warning: package 'foreach' was built under R version 4.2.3
#> Loading required package: iterators
#> Warning: package 'iterators' was built under R version 4.2.3
#> Loading required package: snow
#> Warning: package 'snow' was built under R version 4.2.3
cl <- makeCluster(4, type = "SOCK")
registerDoSNOW(cl)

windowNum <- 250
step <- 125
pt01Frag <- calcAdjFrag(epoch = epoch, window = windowNum, step = step, parallel = TRUE, nSearch=100L,progress = TRUE)

# Fragility result
pt01Frag
#> 
#> Fragility object
#> -----------------------------------------------------------------
#>  Slot       | Class     | DIM/LEN |  Values                     
#> -----------------------------------------------------------------
#>  frag       | matrix    | [84,23] |  0.6310,  0.6618,  0.7484... 
#>  R2         | matrix    | [84,23] |  0.9703,  0.9848,  0.9936... 
#>  lambdas    | numeric   | [23]    |  0.0001,  0.0010,  0.0001... 
#>  startTimes | numeric   | [23]    | -1.0000, -0.8750, -0.7500... 
#>  electrodes | character | [84]    | G1, G2, G3...                
#> ----------------------------------------------------------------- 
#> Use '$attr' to access the data

# Stop the parallel backend
stopCluster(cl)

Predict the SOZ

High value of the fragility correlate with the SOZ. The function estimateSOZ estimates the SOZ from the fragility matrix. By default, it will return the electrode names that are in the top 10% of the fragility values. In this example, pt01Frag$startTimes > 0 corresponds to the time windows that are after the seizure onset. This is not required. However, we found that this usually improves the results.

soz <- estimateSOZ(pt01Frag[ ,pt01Frag$startTimes > 0])
soz
#> [1] "ATT2" "AD2"  "G32"  "ATT1" "AD3"  "G12"  "G13"  "SLT3" "PLT5"

Fragility Statistics

The function fragStat computes quantiles, mean and standard deviation for two electrodes group marked as soz non marked as soz to check that the high value of fragilities. The significant higher distribution of the fragility in the electrode group marked as soz can be used as a prediction of seizure-freedom after resection of these electrodes. The quantile results can be used in a multi-patient study with a machine learning model to detect automatically the SOZ Li et al. (2017) <doi: 10.23919/ACC.2017.7963378>.

stats <- fragStat(pt01Frag, soz)
stats
#> 
#> FragStat object (Summary Statistics by Step)
#> ------------------------------------------------------------
#>  Slot    | Class   | DIM/LEN |  Values                     
#> ------------------------------------------------------------
#>  qmatrix | matrix  | [20,23] |  0.6933,  0.7547,  0.7818... 
#>  meanSOZ | numeric | [23]    |  0.8066,  0.8528,  0.8890... 
#>  meanRef | numeric | [23]    |  0.6738,  0.6780,  0.7161... 
#>  sdSOZ   | numeric | [23]    |  0.1203,  0.1024,  0.0622... 
#>  sdRef   | numeric | [23]    |  0.1730,  0.1849,  0.1476... 
#> ------------------------------------------------------------ 
#> Use '$attr' to access the data

Plot functions for Fragility Matrix

The function plotFragHeatmap produce a heatmap of the fragility and allows to visually check the correlation between the SOZ and sustained high fragility value. Time window around seizure onset [-1:2]s, which includes the following electrodes:

• Electrodes within the surgery region: sozNames
• Electrodes outside of the surgery region: “MLT1”, “MLT2”, “MLT3”, “MLT4”

The parameter sozIndex accepts either electrode names or their indices in the fragility matrix. If sozIndex is provided, the function will display the sozIndex electrodes in blue in the top rows of the heatmap. If sozIndex is not provided, the function will simply generate a heatmap of the fragility matrix.

display <- c(sozNames, "MLT1", "MLT2", "MLT3", "MLT4")
plotFragHeatmap(frag = pt01Frag, sozIndex = sozNames)

plotFragDistribution and plotFragQuantile are similar to the plotFragHeatmap function. They plot the fragility distribution and quantiles, respectively.

plotFragDistribution(frag = pt01Frag, sozIndex = sozNames)

plotFragQuantile(frag = pt01Frag, sozIndex = sozNames)