Uncertainty Quantification¶

Tissue conductivity values are uncertain. Given that, it might be desirable to take this variability into consideration when simulating electric fields.

For that reason, we introduced a way in SimNIBS to easily perform Uncertainty Quantification (UQ) via Generalized Polynomial Chaos Expansion (gPC) for evaluating the effect of uncertainty in conductivity values on the fields of interest.

Before proceeding through this tutorial, please see Saturnino et al., 2018.

Introduction¶

In UQ via gPC, we quantify the uncertainty of the input variables (the conductivities) using a probability distribution. SimNIBS supports two types of conductivity distribution: Normal and Beta.

After defining the uncertainty in the input variables, SimNIBS performs the UQ via gPC. In gPC, we build a polynomial representation of the output variable (such as the electric field), given the input variable. For SimNIBS, we developed an adaptive approach to select the polynomials to obtain fast convergence. At each step, the error is evaluated using a K-means cross-validation scheme. The iteration stops when a tolerance is reached.

For more information on methodology, please see the supplementary material in Saturnino et al., 2018.

Setting-up a Simulation with UQ¶

It is simple to set-up a simulation with UQ. All you need to do is to set-up the distribution_type and distribution_parameters parameters in the COND structures.

Note

As of now, UQ is only supported through scripting. For more information on scripts, please see Scripting Simulations.

Warning

Because conductivities can not be negative, we highly recommend the usage of beta-distributed random variables to represent conductivity uncertainties.

Python¶

''' Example for running a SimNIBS tDCS simulation in Python
    Run with:

    simnibs_python uncertainty_quantification.py

    Copyright (C) 2019 Guilherme B Saturnino
'''
from simnibs import sim_struct, run_simnibs

S = sim_struct.SESSION()
S.subpath = 'm2m_ernie'
S.pathfem = 'tdcs_uq'

tdcs = S.add_tdcslist()
tdcs.currents = [0.001, -0.001]

# Set-up the electrodes
mc_electrode = tdcs.add_electrode()
mc_electrode.channelnr = 1
mc_electrode.centre = 'C3'
mc_electrode.shape = 'rect'
mc_electrode.dimensions = [50, 50]
mc_electrode.thickness = 4

so_electrode = tdcs.add_electrode()
so_electrode.channelnr = 2
so_electrode.centre = 'AF4'
so_electrode.shape = 'rect'
so_electrode.dimensions = [50, 70]
so_electrode.thickness = 4

# Set-up the uncertain conductivities
# White Matter
tdcs.cond[0].distribution_type = 'beta'
tdcs.cond[0].distribution_parameters = [3, 3, .1, .4]
# Gray Matter
tdcs.cond[1].distribution_type = 'beta'
tdcs.cond[1].distribution_parameters = [3, 3, .1, .6]
# Compact Bone
tdcs.cond[6].distribution_type = 'beta'
tdcs.cond[6].distribution_parameters = [3, 3, 0.001, 0.012]

run_simnibs(S)

MATLAB¶

% example script that runs an Uncertainty Quantification analysis
% 
% G. Saturnino, 2019
%

%% General information

S = sim_struct('SESSION');
S.subpath = 'm2m_ernie'; % subject folder
S.pathfem = 'tdcs_uq'; %Folder for the simulation output

%% Define tDCS simulation
S.poslist{1} = sim_struct('TDCSLIST');
S.poslist{1}.currents = [0.001, -0.001]; % Current flow though each channel (mA)

%First Electrode
S.poslist{1}.electrode(1).channelnr = 1; % Connect the electrode to the first channel
S.poslist{1}.electrode(1).centre = 'C3'; % Place it over C3
S.poslist{1}.electrode(1).shape = 'rect'; %Rectangular electrode
S.poslist{1}.electrode(1).dimensions = [50, 50]; % Dimension in mm
S.poslist{1}.electrode(1).thickness = 4; % 4 mm thickness

%Second Electrode
S.poslist{1}.electrode(2).channelnr = 2;
S.poslist{1}.electrode(2).centre = 'AF4';
S.poslist{1}.electrode(2).shape = 'rect';
S.poslist{1}.electrode(2).dimensions = [50, 70];
S.poslist{1}.electrode(2).thickness = 4;

% Set-up the uncertain conductivities
% White Matter
S.poslist{1}.cond(1).distribution_type = 'beta';
S.poslist{1}.cond(1).distribution_parameters = [3, 3, .1, .4];
% Gray Matter
S.poslist{1}.cond(2).distribution_type = 'beta';
S.poslist{1}.cond(2).distribution_parameters = [3, 3, .1, .6];
% Compact Bone
S.poslist{1}.cond(7).distribution_type = 'beta'
S.poslist{1}.cond(7).distribution_parameters = [3, 3, 0.001, 0.012];

%% Run Simulation
run_simnibs(S)

Output Files¶

In the output folder of the examples above (the tdcs_uq/ folder) we have the output files

ernie_TDCS_1_electrodes.msh
Head mesh (Gmsh format) with the electrodes
ernie_TDCS_1_gpc.msh
Uncertainty quantification result in Gmsh format. There we have the data sets:
- {Field Name}_mean
  Mean of the probability distribution describing the uncertainty of the given field
- {Field Name}_std
  Standard Deviation of the probability distribution describing the uncertainty of the given field
- {Field Name}_sensitivity_{N}
  Derivative-based sensitivity for the tissue number N. See Saturnino et al., 2018
- {Field Name}_sobol_{N}
  Sobol coefficient for the tissue number N. See Saturnino et al., 2018
ernie_TDCS_1_gpc.hdf5: HDF5 file containing partial information on the gPC expansion and partial results. The main datasets there are:
- gpc_object
  Information to reconstruct the gPC regressions
- cond
  Conductivity information
- mesh
  Original mesh (or with electrodes for tDCS simulations). Has the mesh/data_matrices/v_samples data set with the electric potential values at all mesh nodes for all iterations of the adaptive algorithm.
- mesh_roi
  Mesh cropped to the ROI. Has the mesh/data_matrices/Field Name_samples data sets with the field values (at nodes for v or elements for other quantities) for all iterations of the adaptive algorithm.
Note

meshes stored in HDF5 files can be read with the simnibs.msh.Msh.read_hdf5() class or with the mesh_load_hdf5 function in MATLAB

More Options¶

In Python, it is also possible to call lower-level functions to set more options for the UQ run and the post-processing of results.

In the example below, we set-up a UQ TMS problem with the ROI being the whole brain (tissues 1 and 2) and with a tolerance of 0.1.

import os

import simnibs
from simnibs.simulation import gpc

## Define the TMS simulation
tms = simnibs.sim_struct.TMSLIST()
tms.fnamecoil = os.path.join('legacy_and_other','Magstim_70mm_Fig8.ccd')
tms.mesh = simnibs.read_msh(os.path.join('m2m_ernie','ernie.msh'))

# We need to set the EEG cap as we are giving the coil
# Coordinates as EEG positions
sub_files = simnibs.SubjectFiles(subpath='m2m_ernie')
tms.eeg_cap = sub_files.eeg_cap_1010

# Define the coil position
pos = tms.add_position()
pos.centre = 'C3'
pos.pos_ydir = 'CP3'
pos.distance = 4

# Define the uncertain conductivities
tms.cond[0].distribution_type = 'beta'
tms.cond[0].distribution_parameters = [3, 3, .1, .4]
tms.cond[1].distribution_type = 'beta'
tms.cond[1].distribution_parameters = [3, 3, .1, .6]
tms.cond[6].distribution_type = 'beta'
tms.cond[6].distribution_parameters = [3, 3, 0.003, 0.012]

# Run the UQ calling with White and Gray matter as an ROI and tolerance of 1e-2
gpc.run_tms_gpc(tms, 'tms_gpc/TMS_UQ', tissues=[1, 2], eps=1e-2)

Secondary Quantities¶

It is also possible to calculate secondary quantities, such as the 99th percentile of the electric field magnitude.

import numpy as np
import simnibs
from simnibs.simulation import gpc


fn_hdf5 = 'tdcs_uq/ernie_TDCS_1_gpc.hdf5'
# Read the regression object from the HDF5 file
regression = gpc.gPC_regression.read_hdf5(fn_hdf5)
# Read the mesh ROI from the HDF5 file
mesh = simnibs.Msh.read_hdf5(fn_hdf5, 'mesh_roi')
# Define the function to be used for the expansion

def percentile_99(Es):
    # The function will receive the electric field in a format
    # N_sims x N_elm x 3
    # for each simulation, we calculate the 99th percentile
    prc = np.zeros(Es.shape[0])
    for i, E in enumerate(Es):
        prc[i] = simnibs.ElementData(E, mesh=mesh).get_percentiles(99)[0]

    return prc

# Calculate the gPC coefficients
gpc_coeffs = regression.expand_quantity(percentile_99)

print("99th Percentile")
print("Mean Value: ", regression.mean(gpc_coeffs))
print("Standard Deviation: ", regression.std(gpc_coeffs))

# Draw 1000 samples for the 99th percentile
samples = regression.MC_sampling(gpc_coeffs, 1000)[1]

Acknowledgements¶

We would like to thank Konstantin Weise and Thomas Knoesche for the support in implementing the gPC in SimNIBS and supplying us with an early version of the pygpc library.

Uncertainty Quantification¶

Introduction¶

Setting-up a Simulation with UQ¶

Python¶

MATLAB¶

Output Files¶

More Options¶

Secondary Quantities¶

Acknowledgements¶

Further Reading¶

Navigation

Related Topics