The hardware and bandwidth for this mirror is donated by METANET, the Webhosting and Full Service-Cloud Provider.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]metanet.ch.

manuscript_reproducible

Anne-Cecile Lesage, PhD, Oliver Zhou, BA, Jiefei Wang, PhD

2025-12-26

This vignette shows how to reproduce results from the companion software paper of the EZFragility CRAN R package for the first seizure of patients 01 and 26. This allows to check how well the EZFragility package reproduces results from the original neural fragility paper

Load the package and data

library(EZFragility)
library(Epoch)
library(ggplot2)
library(ggtext)
library(gsignal)
#> 
#> Attaching package: 'gsignal'
#> The following object is masked from 'package:Epoch':
#> 
#>     resample
#> The following objects are masked from 'package:stats':
#> 
#>     filter, gaussian, poly

Data preprocess (Not included in the package)

To allow fast and easy exploration of fragility results, we preprocessed the multipatient data from the OpenNeuro Fragility Data Set with RAVE 2.0 . Main steps included Notch filtering to remove power line interference, re-referencing, and epoching based on seizure onset time to isolate relevant segments of [-30:30s].

Download Data

To access the preprocessed data, you can use the EpochDownloader class from the Epoch package. The following code downloads the data and lists the available datasets.

dl<-EpochDownloader(progress = FALSE)
names(dl)
#>  [1] "FragilityData_subpt2_3"   "FragilityData_subpt13_2" 
#>  [3] "FragilityData_subpt15_1"  "FragilityData_subjh105_2"
#>  [5] "FragilityData_subpt12_1"  "FragilityData_subjh103_1"
#>  [7] "FragilityData_subpt7_2"   "FragilityData_subpt11_2" 
#>  [9] "FragilityData_subpt13_1"  "FragilityData_subpt8_2"  
#> [11] "FragilityData_subpt2_2"   "FragilityData_subjh103_2"
#> [13] "FragilityData_subjh105_1" "FragilityData_subpt6_3"  
#> [15] "FragilityData_subpt6_2"   "FragilityData_subpt8_1"  
#> [17] "FragilityData_subjh101_2" "FragilityData_subjh103_3"
#> [19] "FragilityData_subpt16_1"  "FragilityData_subpt01_1" 
#> [21] "FragilityData_subpt01_4"  "FragilityData_subpt12_2" 
#> [23] "FragilityData_subpt3_1"   "FragilityData_subpt10_1" 
#> [25] "FragilityData_subpt11_4"  "FragilityData_subpt2_1"  
#> [27] "FragilityData_subpt11_3"  "FragilityData_subpt11_1" 
#> [29] "FragilityData_subpt14_2"  "FragilityData_subpt15_2" 
#> [31] "FragilityData_subpt10_3"  "FragilityData_subpt3_2"  
#> [33] "FragilityData_subpt01_3"  "FragilityData_subjh101_1"
#> [35] "FragilityData_subjh105_3" "FragilityData_subpt10_2" 
#> [37] "FragilityData_subpt15_4"  "FragilityData_subpt8_3"  
#> [39] "FragilityData_subpt01_2"  "FragilityData_subpt6_1"  
#> [41] "FragilityData_subpt15_3"  "FragilityData_subpt16_3" 
#> [43] "FragilityData_subjh105_4" "FragilityData_subpt14_1" 
#> [45] "FragilityData_subpt16_2"  "FragilityData_subpt7_3"  
#> [47] "FragilityData_subpt17_2"

The preprocessed voltage data from patient pt01 seizure 1 and pt26 seizure 1 can be loaded by

pt01sz1 <- dl$FragilityData_subpt01_1
pt01sz1
#> Epoch Object @ 1000 Hz: 
#>  Time -30       -29.999   -29.998   -29.997   -29.996   -29.995   -29.994   
#>  G1   -340001.1 -337468.1 -360006   -371602.5 -381722.9 -388522.7 -409074   ... 
#>  G2   -258822.7 -260679   -262803.7 -258883.7 -254481.7 -252272.8 -258038.6 ... 
#>  G3   -239279.3 -233117.6 -223820.4 -206524.9 -194307.3 -175271.1 -163560.5 ... 
#>  G4   123333.6  116130.2  121380.4  133184.7  156112.5  191592.9  235748.3  ... 
#>  G7   63804.24  59270.13  51457.29  36245.07  17488.84  4907.711  -14826.07 ... 
#>  G8   145233.7  132424.9  126525    115485.8  112664.6  106142.3  95139.39  ... 
#>  G9   -245320.8 -238745.9 -241513.8 -253517.7 -264246.5 -271382   -270364.3 ... 
#>  G10  -495154.8 -480099.8 -473566.2 -468136   -469964.3 -477490.8 -485846.7 ... 
#>  G13  209178.2  211959.4  214446.7  230023.3  251683.6  261378.3  274732.9  ... 
#>  G14  170761.3  165173.8  156977.7  156472.7  160976.6  162013.1  158694.1  ... 
#> ...
#> [82 rows x 60001 cols]
#> rowData: [3 vars] name, soz, resected 
#> metaData: [14 vars] patient, type, timeWindow, seizureTime, totalRuns, run, ... 
#> Use tblData, rowData, colData, metaData to get the data

pt26sz1 <- dl$FragilityData_subjh103_1
pt26sz1
#> Epoch Object @ 1000 Hz: 
#>  Time -30       -29.999  -29.998   -29.997   -29.996   -29.995   -29.994   
#>  ABT1 1201953   1185980  1173963   1165929   1161175   1154228   1145490   ... 
#>  ABT2 420150    415039.8 407836.8  399440    392928.7  384582.9  375171.1  ... 
#>  ABT3 362938.1  361439.3 358149.9  351630.6  346282.1  339042.7  331386.2  ... 
#>  ABT4 -98432.64 -96252.9 -96281.68 -98979.26 -98652.86 -104023.3 -114738.3 ... 
#>  MBT1 378468.9  377128.7 375202.9  374722.7  376949.6  376394.8  377590.7  ... 
#>  MBT2 517111.9  508906.5 499356.3  491326    484962.9  477177.7  468958.4  ... 
#>  MBT3 500315.6  496029.1 489656.4  481150    474300.6  463975.6  450947.4  ... 
#>  MBT4 79912.36  75784.27 69818.25  65409     65214.41  67384.24  67453.96  ... 
#>  PBT3 612871.8  605579.9 596700.3  587817.1  579319.7  568749.6  557960.1  ... 
#>  PBT4 267135.9  261987.7 258530.2  255491.5  256614.3  256436.3  256095.7  ... 
#> ...
#> [87 rows x 60001 cols]
#> rowData: [3 vars] name, soz, resected 
#> metaData: [14 vars] patient, type, timeWindow, seizureTime, totalRuns, run, ... 
#> Use tblData, rowData, colData, metaData to get the data

Remove Artifacts

The following function applies a band-pass filter voltage between 0.5 and 150 frequency with fourth-order Butterworth filter to the Epoch objects to remove high-frequency artifacts.

butterworthFilter <- function(epoch, lowpass=0.5, highpass=150) {
  order <- 4

  sampling_freq <- metaData(epoch)$samplingRate
  nyquist_freq <- sampling_freq / 2

  normalized_freqs <- c(lowpass, highpass) / nyquist_freq
  
  filter_type <- "pass"
  butter_filter <- gsignal::butter(
    n = order,
    w = normalized_freqs,
    type = filter_type)

  # Apply filter to epoch data
  mat <- tblData(epoch)

  # Apply zero-phase filtering (filtfilt) to each row
  filtered_data <- gsignal::filtfilt(
    filt = butter_filter,
    x = t(mat))

  filtered_data <- t(filtered_data)
  tblData(epoch) <- filtered_data

  epoch
}

We apply the filter to the epochs and crop the data to the relevant time window of [-10:10s] around seizure onset:

pt01sz1Cropped <- pt01sz1 |>
    crop(start=-10, end=10) |>
    butterworthFilter()
pt01sz1Cropped
#> Epoch Object @ 1000 Hz: 
#>  Time -10       -9.999    -9.998    -9.997    -9.996    -9.995    -9.994    
#>  G1   -97900.28 -92900.03 -90766.22 -92447.49 -96380.57 -99523.05 -99226.33 ... 
#>  G2   -56052.66 -40263.26 -29888.38 -27397.71 -31305.75 -37609.67 -42506.36 ... 
#>  G3   -15517.31 -17632.76 -18318.98 -16737.49 -13053.21 -8498.531 -4982.956 ... 
#>  G4   -49128.8  -50584.73 -53861.98 -60145.2  -69759.44 -82404.06 -97606.44 ... 
#>  G7   -37898.06 -31519.74 -26516.46 -23804.4  -23456.72 -24600.91 -25770.64 ... 
#>  G8   52495.89  39420.95  29423.77  24669.29  26019.12  33303.49  45919.68  ... 
#>  G9   -18412.16 -34190.81 -49819.09 -63998.55 -73961.42 -76265.79 -68360.34 ... 
#>  G10  -51285.93 -49613.51 -51586.78 -58840.26 -70102.13 -81806.09 -89814.83 ... 
#>  G13  69408.1   64320.8   57154.87  47208.05  35591.58  24488.22  15864.63  ... 
#>  G14  146480.8  135457.6  127275    123597.7  124438.5  128573.1  134397.6  ... 
#> ...
#> [82 rows x 20001 cols]
#> rowData: [3 vars] name, soz, resected 
#> metaData: [14 vars] patient, type, timeWindow, seizureTime, totalRuns, run, ... 
#> Use tblData, rowData, colData, metaData to get the data

pt26sz1Cropped <- pt26sz1 |>
    crop(start=-10, end=10) |>
    butterworthFilter()
pt26sz1Cropped
#> Epoch Object @ 1000 Hz: 
#>  Time -10       -9.999    -9.998    -9.997    -9.996    -9.995    -9.994    
#>  ABT1 27127.52  29593.17  32935.09  36969.52  40467.82  42209.74  42090.3   ... 
#>  ABT2 -39902.05 -42339.29 -44080.98 -44875.34 -45038.9  -45237.56 -46132.1  ... 
#>  ABT3 -43984.3  -50470.35 -57409.72 -64842.17 -72191.7  -78422.62 -82524.66 ... 
#>  ABT4 -2441.398 -4379.817 -5206.001 -4434.294 -2394.185 74.10756  2122.067  ... 
#>  MBT1 -19055.05 -23234.11 -26628.32 -29037.06 -30875.35 -32640.04 -34349.15 ... 
#>  MBT2 -12796.88 -14464.27 -15915.51 -17022.43 -17756.34 -18125.28 -18110.47 ... 
#>  MBT3 1295.278  -1771.742 -4850.438 -7768.158 -10092.93 -11267.79 -10935.69 ... 
#>  MBT4 -44316.8  -46028.58 -47239.21 -47803.94 -47914.95 -47781.69 -47423.46 ... 
#>  PBT3 -12199.13 -10197.99 -7971.241 -5618.192 -3536.198 -2119.899 -1527.249 ... 
#>  PBT4 -2833.469 899.8083  4697.105  8171.53   10481.06  10777.59  8722.417  ... 
#> ...
#> [87 rows x 20001 cols]
#> rowData: [3 vars] name, soz, resected 
#> metaData: [14 vars] patient, type, timeWindow, seizureTime, totalRuns, run, ... 
#> Use tblData, rowData, colData, metaData to get the data

Visualize Voltage Plot (Figure 2.b)

We show the voltage plot using plot(epoch) for a selected subset of electrodes for patient 01 and 26. We will highlight the SOZ electrodes in red and non-SOZ electrodes in black.

# A helper function to visualize SOZ electrodes in red and non-SOZ in black
visualSOZ <- function(epoch, sozNames) {
  p <- plot(epoch, gap = 4, timeResolution = 512)

  elecColor <- rep("black", nrow(epoch))
  elecColor[rownames(epoch) %in% sozNames] <- 'red'
  elecColor <- rev(elecColor) # match the electrode order in the plot

  p +
    geom_vline(xintercept = 0, color = "black", linetype = "dashed", linewidth = 1)+
    theme(axis.text.y = element_markdown(colour = elecColor))
}

To visualize patient 01 seizure 1:

pt01Subset <- c(
    "ATT1", "ATT2", "AD1", "AD2", "AD3", 
    "AD4", "PD1", "PD2", "PD3", "PD4", 
    "MLT1", "MLT2", "MLT3", "MLT4")
pt01sozName <- c(
    "ATT1", "ATT2", "AD1", "AD2", "AD3", 
    "AD4", "PD1", "PD2", "PD3", "PD4"
)
pt01sz1Reordered <- pt01sz1Cropped[pt01Subset, ]
visualSOZ(pt01sz1Reordered, pt01sozName)

To visualize patient 26 seizure 1:

pt26Subset <- c(
    "ABT1", "ABT2", "RAD1", "RAD2", "RAD3", "RAD4", "RAD5", "RAD6", 
    "RAD7", "RHD1", "RHD2", "RHD3", "RHD4", "RHD5", "RHD6", "RHD7",
    "RHD8", "RHD9", "RTG29", "RTG30", "RTG31", "RTG32", "RTG40",
    "RTG48")
pt26sozName <- c(
    "RAD1", "RAD2", "RAD3", "RAD4", "RAD5", "RAD6", "RAD7", "RHD1", 
    "RHD2", "RHD3", "RHD4", "RHD5", "RHD6", "RHD7", "RHD8", "RHD9",
    "RTG40", "RTG48")
pt26sz1Reordered <- pt26sz1Cropped[pt26Subset, ]
visualSOZ(pt26sz1Reordered, pt26sozName)

Compute the Fragility Matrix

The following code computes the fragility matrix using all electrodes and store the results in the Fragility class object pt01sz1Frag

library(doSNOW)
# compute fragility

cl <- makeCluster(parallel::detectCores(), type = "SOCK")
registerDoSNOW(cl)

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

pt26sz1Frag <- calcAdjFrag(epoch = pt26sz1Cropped, window = windowNum, step = step, parallel = TRUE, nSearch=100L, progress = FALSE)

# Stop the parallel backend
stopCluster(cl)

Fragility Heatmap (Figure 2.a)

We plot the fragility heatmap using plot(frag) with the same display options as the previous voltage plot. Looking at both plots allows to check correlation between soz patterns

# A helper function to visualize fragility heatmap with SOZ electrodes in red and non-SOZ in black
fragHeatmap <- function(frag, sozNames, ranked=FALSE) {
  startTimes <- frag$startTimes

  indexsz <- which(abs(startTimes)<=0.01)
  elecColor <- rep("black", length(frag$electrodes))
  elecColor[frag$electrodes%in% sozNames] <- 'red'
  elecColor <- rev(elecColor)

  plot(frag, ranked=ranked) +
    geom_vline(xintercept = indexsz, color = "black", linetype = "dashed", linewidth = 1) +
    theme(
      axis.text.y = element_markdown(colour = elecColor)
    )
}

For patient 01, we can see that the SOZ electrodes identified by clinicians have high fragility values compared to the non-SOZ electrodes. The regions corresponding to the SOZ electrodes were resected during surgery, and patient 01 was a surgical success.

pt01sz1FragReordered <- pt01sz1Frag[pt01Subset]
fragHeatmap(pt01sz1FragReordered, pt01sozName)

For patient 26, many non-SOZ electrodes also had high fragility values. The regions corresponding to those high-fragility non-SOZ electrodes were not resected during surgery, and patient 26 was a surgical failure.

pt26sz1FragReordered <- pt26sz1Frag[pt26Subset]
fragHeatmap(pt26sz1FragReordered, pt26sozName)

Mean and standard deviation statistics for SOZ and non SOZ electrode group (Figure 3.a)

We use plotFragDistribution(frag) function from the EZFragility package to visualize the mean and standard deviation statistics for the SOZ and non-SOZ electrode groups. This plot shows that the fragility biomarker statistics are respectively significantly/not significantly higher in the soz labeled group correlated with the ground truth success/failure outcome for patient 01 and patient 26.

# A helper function to plot fragility distribution for SOZ and non-SOZ groups
fragDist <- function(frag, sozNames) {
    timeWindows <- frag$startTimes
    timeIdx <- which(timeWindows >= -5 & timeWindows <= 10)
    frag <- frag[, timeIdx]
    plotFragDistribution(frag = frag, groupIndex = sozNames, bandType="SEM", rollingWindow = 1) +
        geom_vline(xintercept = 0, color = "black", linetype = "dashed", linewidth = 1)
}

Patient 01 seizure 1 (p 1481 Extended Fig.4)

fragDist(pt01sz1Frag[pt01Subset], pt01sozName)

Patient 26 seizure 1 (p 1481 Extended Fig.4)

fragDist(pt26sz1Frag[pt26Subset], pt26sozName)

Distribution of Fragility Values (Figure 3.b)

We use plotFragQuantile(frag) function from the EZFragility package to visualize the distribution of fragility values across all time windows for the SOZ and non-SOZ electrode groups. For patient 01:

plotFragQuantile(pt01sz1Frag[pt01Subset])

For patient 26:

plotFragQuantile(pt26sz1Frag[pt26Subset])

Comparing different lambda values (Figure 4)

Since the fragility computation relies on L2-norm regularization, we can explore how different lambda values impact the fragility results. We compute fragility matrices for three different lambda values (1e-4, 1e-3, 1e-2) and visualize the fragility heatmaps for patient 01 seizure 1.

library(doSNOW)
cl <- makeCluster(parallel::detectCores(), type = "SOCK")
registerDoSNOW(cl)

windowNum <- 250
step <- 125

pt01sz1Frag_lambda1 <- calcAdjFrag(epoch = pt01sz1Cropped, window = windowNum, step = step, parallel = TRUE, nSearch=100L, progress = FALSE, lambda=1e-4)


pt01sz1Frag_lambda2 <- calcAdjFrag(epoch = pt01sz1Cropped, window = windowNum, step = step, parallel = TRUE, nSearch=100L, progress = FALSE, lambda=1e-3)

pt01sz1Frag_lambda3 <- calcAdjFrag(epoch = pt01sz1Cropped, window = windowNum, step = step, parallel = TRUE, nSearch=100L, progress = FALSE, lambda=1e-2)


# Stop the parallel backend
stopCluster(cl)

We can see that the choice of lambda impacts the fragility values, but the overall patterns seems relatively stable. However, it should be noted that there is no mathmatical guarantee that the fragility values are valid when the connection matrix is unstable.

fragHeatmap(pt01sz1FragReordered, pt01sozName, ranked=FALSE) + ggtitle("λ = Auto-selected")

fragHeatmap(pt01sz1Frag_lambda1[pt01Subset], pt01sozName, ranked=FALSE) + ggtitle("λ = 1e-4")

fragHeatmap(pt01sz1Frag_lambda2[pt01Subset], pt01sozName, ranked=FALSE) + ggtitle("λ = 1e-3")

fragHeatmap(pt01sz1Frag_lambda3[pt01Subset], pt01sozName, ranked=FALSE) + ggtitle("λ = 1e-2")

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.