Neural model state and parameter estimation from data
Copyright (C) 2018 Dean R. Freestone, Philippa J. Karoly - All Rights Reserved
You may use, distribute and modify this code under the terms of the GNU General Public license v2.0. You can find a copy of this license in file LICENSE.md. Note than when distributing derived works, the source code of the work must be made available under the same license.
Figures in this repository are licensed under Creative Commons with conditions:
- Attribution: please cite Seizure Pathways: A model-based investigation, P.J. Karoly, L. Kuhlmann, D. Soudry, D.B. Grayden, M.J. Cook and D.R. Freestone (2018) PLoS Computational Biology (https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1006403)
To cite this code in academic research please reference: Seizure Pathways: A model-based investigation, P.J. Karoly, L. Kuhlmann, D. Soudry, D.B. Grayden, M.J. Cook and D.R. Freestone (2018) PLoS Computational Biology (https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1006403)
To cite the data in academic researh please reference: Melbourne University NeuroVista Seizure Prediction Data (https://doi.org/10.26188/5b6a999fa2316) and Cook, M. J., O'Brien, T. J., Berkovic, S. F., Murphy, M., Morokoff, A., Fabinyi, G., ... & Hosking, S. (2013). Prediction of seizure likelihood with a long-term, implanted seizure advisory system in patients with drug-resistant epilepsy: a first-in-man study. The Lancet Neurology, 12(6), 563-571.
Source code is located in src folder
- set_params.m: A function that sets the parameters of the neural mass model (Jansen & Rit 1995)
- g.m: A function that computes the erf sigmoid
- prop_mean_covariance.m: A function that computes the posterior mean and covariance of the state/parameter distribution propagated through the neural mass model
Example code is located in Example folder
- Example.m: Runs the estimation for a single seizure
- generateData.m: Runs the model at different values of input and generates simulated data
One example seizure is provided in data folder. Additional seizure data can be downloaded online from the Epilepsy Ecosystem. Data are available as “Melbourne Seizure Prediction Trial Seizure Data” from https://www.epilepsyecosystem.org/howitworks/#data. After registering you will be emailed an invitation to download the data from figshare. The data are licensed under a Creative Commons license with conditions on ATTRIBUTION, NON-COMMERCIAL use and SHARE-ALIKE. Users will be required to create an account and sign a terms of use agreement that requires no commercial use, and restricts all works derived from the data to be made publicly available under a Creative Commons license.
Data is a single .mat file with a variable, Seizure that has dimension T x N, where T is the number of samples, and N is the number of electrode channels (16). The seizure onset is 5 minutes from the start of the data and seizure offset is 1 minute from the end of the data. Data is sampled at 400Hz.
These figures relate specifically to the results presented in Karoly et al (2018). We provide additional figures for connectivity parameter estimation and signal energy showing all 16 channels.
If you are not familiar with Kalman filtering, a review of the recommended resources (or similar) is strongly advised before implementing this code. Density filters can be plagued by numerical instability and in practise fine-tuning of filter parameters is often required to run the filter on real-world data. We provide a non-exhaustive list of some gotchas and heuristics that we have observed over the years.
- Data range is inconsistent with the model (scale your data so it lies not to far from the bounds of what your model can simulate)
- State/parameter values differ by many orders of magnitude (scale your model)
- Mismatched DC between model and data (add one constant offset parameter to the estimate)
- Covariance, P becomes asymmetric (enforce symmetry with P = (P + P') / 2)
- Covariance becomes too small, breaking the numerics of the filter (add a very small amount to the diagonal of Q, initialise P larger)
- Estimates don't converge (try using annealing to gradually increase R)
- Too many parameters to estimate (maybe you need to lump them together)
More information can be found in the following refrences
-
Kuhlmann, L., Freestone, D. R., Manton, J. H., Heyse, B., Vereecke, H. E., Lipping, T., ... & Liley, D. T. (2016). Neural mass model-based tracking of anesthetic brain states. NeuroImage, 133, 438-456.
-
Optimal State Estimation: Kalman, H Infinity, and Nonlinear Approaches, Dan Simon and the authors example code!
-
Neural Control Engineering: The Emerging Intersection Between Control Theory, Steve Schiff (esp Chapter 2 and Chapter 5) and all the example code