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 in SimNIBS a way 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 scipts, 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.fnamehead = 'ernie.msh'
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]
# Skull
tdcs.cond[3].distribution_type = 'beta'
tdcs.cond[3].distribution_parameters = [3, 3, 0.003, 0.012]

run_simnibs(S)

MATLAB

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

%% General information

S = sim_struct('SESSION');
S.fnamehead = 'ernie.msh'; % head mesh
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];
% Skull
S.poslist{1}.cond(4).distribution_type = 'beta';
S.poslist{1}.cond(4).distribution_parameters = [3, 3, 0.003, 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
from simnibs import sim_struct, msh, file_finder
from simnibs.simulation import gpc

## Define the TMS simulation
tms = sim_struct.TMSLIST()
tms.fnamecoil = 'Magstim_70mm_Fig8.nii.gz'
tms.mesh = msh.read_msh('ernie.msh')

# We need to set the EEG cap as we are giving the coil
# Coordinates as EEG positions
sub_files = file_finder.SubjectFiles('ernie.msh')
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[3].distribution_type = 'beta'
tms.cond[3].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 norm

import numpy as np
from simnibs import msh
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 = msh.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] = msh.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]
breakpoint()

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.