TDCS Optimization

SimNIBS 3 is able of optimizing TDCS montages for maximal focality. There are two steps for performing the optimizations.

  1. Create a leadfield.

  2. Set-up an optimization based on the leadfield.

Leadfield Calculations

The leadfield is calculated by first placing many electrodes in the head accordingly with an EEG cap, for example, and afterwards calculating the field caused by each electrode individually, keeping a constant return electrode.

With this simulations, we form a matrix which can then be used to calculate the electric field caused by any combination of electrodes. As the leadfield involves calculating many electric fields, it takes some time to run, but usually no more than 1h.

To calculate leadfields, the user must write a small Python of Matlab script using the TDCSLEADFIELD structure. Unless specified differently, the leadfields are calculated for EEG10-10 electrode positions.

  • Python

    ''' Example of a SimNIBS tDCS leadfield in Python
        Run with:
    
        simnibs_python leadfield.py
    
        Copyright (C) 2019 Guilherme B Saturnino
        
        place script in the “ernie” folder of the example dataset
    '''
    from simnibs import sim_struct, run_simnibs
    
    tdcs_lf = sim_struct.TDCSLEADFIELD()
    # head mesh
    tdcs_lf.fnamehead = 'ernie.msh'
    # output directory
    tdcs_lf.pathfem = 'leadfield'
    
    
    # Uncoment to use the pardiso solver
    #tdcs_lf.solver_options = 'pardiso'
    # This solver is faster than the default. However, it requires much more
    # memory (~12 GB)
    
    
    run_simnibs(tdcs_lf)
    

  • MATLAB

    % Example of a SimNIBS tDCS leadfield 
    % Copyright (C) 2019 Guilherme B Saturnino
    
    % place script in the 'ernie' folder of the example dataset
    
    tdcs_lf = sim_struct('TDCSLEADFIELD');
    % Head mesh
    tdcs_lf.fnamehead = 'ernie.msh';
    % Output directory
    tdcs_lf.pathfem = 'leadfield';
    
    % Uncomment to use the pardiso solver
    %tdcs_lf.solver_options = 'pardiso';
    % This solver is much faster than the default. However, it requires much more
    % memory (~12 GB)
    
    run_simnibs(tdcs_lf)
    

Here, we are calculating a leadfield in the ernie example data set, based on the EEG 10-10 system.

Note

For more options in running tDCS leadfields, please see the documentation at TDCSLEADFIELD.

This script will create the files in the leadfield folder:

  • ernie_electrodes_EEG10-10_UI_Jurak_2007.msh: File with the head mesh and electrodes placed on top. Can be used for checking electrode positions. Beware that some electrodes will be of the same color as the skin.

  • ernie_ROI.msh: The region of interest where the electric field values were recorded. In the example, the middle gray matter surface.

  • ernie_leadfield_EEG10-10_UI_Jurak_2007.hdf5: HDF5 file with the leadfield matrix, orgenized with the datasets:

    • mesh_electrodes: Full mesh with the electrodes

    • mesh_leadfield: Mesh with the region of interest only (in the example, the middle gray matter surface). Has the mesh_leadfield/leadfields/tdcs_leadfield data set with the leadfield. In the attributes, you can see information such as the reference electrode, the field (E or J), units, currents and so on.

    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

Optimization

Now, we will use the leadfield to optimize the electric field at a given target.

Simple Optimiztion

The first step is to set TDCSoptimize structure. In this structure, we need to select the leadfield we will use for the optimization, a name for the optimization problem, safety constraints and limit the number of electrodes.

Afterwards, we need to define an optimization target using an TDCStarget structure.

  • Python

    ''' Example of a SimNIBS tDCS optimization in Python
        Run with:
    
        simnibs_python tdcs_optimize.py
    
        Copyright (C) 2019 Guilherme B Saturnino
    '''
    
    
    import simnibs
    
    # Initialize structure
    opt = simnibs.opt_struct.TDCSoptimize()
    # Select the leadfield file
    opt.leadfield_hdf = 'leadfield/ernie_leadfield_EEG10-10_UI_Jurak_2007.hdf5'
    # Select a name for the optimization
    opt.name = 'optimization/single_target'
    
    # Select a maximum total current (in A)
    opt.max_total_current = 2e-3
    # Select a maximum current at each electrodes (in A)
    opt.max_individual_current = 1e-3
    # Select a maximum number of active electrodes (optional)
    opt.max_active_electrodes = 8
    
    # Define optimization target
    target = opt.add_target()
    # Position of target, in subject space!
    # please see tdcs_optimize_mni.py for how to use MNI coordinates
    target.positions = [-55.4, -20.7, 73.4]
    # Intensity of the electric field (in V/m)
    target.intensity = 0.2
    
    # Run optimization
    simnibs.run_simnibs(opt)
    

  • MATLAB

    % Example of a SimNIBS tDCS optimization in MATLAB
    %
    % Copyright (C) 2019 Guilherme B Saturnino
    
    
    % Initialize structure
    opt = opt_struct('TDCSoptimize');
    % Select the leadfield file
    opt.leadfield_hdf = 'leadfield/ernie_leadfield_EEG10-10_UI_Jurak_2007.hdf5';
    % Select a name for the optimization
    opt.name = 'optimization/single_target';
    
    % Select a maximum total current (in A)
    opt.max_total_current = 2e-3;
    % Select a maximum current at each electrodes (in A)
    opt.max_individual_current = 1e-3;
    % Select a maximum number of active electrodes (optional)
    opt.max_active_electrodes = 8;
    
    % Define optimization target
    % Position of target, in subject space!
    % please see tdcs_optimize_mni.m for how to use MNI coordinates
    opt.target.positions = [-55.4, -20.7, 73.4];
    % Intensity of the electric field (in V/m)
    opt.target.intensity = 0.2;
    % Run optimization
    run_simnibs(opt);
    

Note

For more information see the documentation for TDCSoptimize and TDCStarget.

Output files

The optimization outputs:

  • name.csv: comma separated values (CSV) files with optimal current values at each electrode (in A)

  • name_electrodes.geo: Gmsh .geo file for visualizing electrodes and currents

  • name.msh: Gmsh .msh file with the target and the optimized electric field in the ROI.

  • name_summary.txt: Some summary quantities about the optimization

Maximizing intensity

To maximize intensity at the target, disregarding field focality, simply use a large value for the target intensity.

  • Python

    ''' Example of optimizaion maximizing intensity in the target
    
        Copyright (C) 2019 Guilherme B Saturnino
    '''
    import simnibs
    
    opt = simnibs.opt_struct.TDCSoptimize()
    opt.leadfield_hdf = 'leadfield/ernie_leadfield_EEG10-10_UI_Jurak_2007.hdf5'
    opt.name = 'optimization/single_target'
    
    opt.max_total_current = 2e-3
    opt.max_individual_current = 1e-3
    opt.max_active_electrodes = 8
    
    target = opt.add_target()
    target.positions = [-55.4, -20.7, 73.4]
    # Set the intensity to a large value (e.g, 100)
    target.intensity = 100
    
    # Run optimization
    simnibs.run_simnibs(opt)
    

  • MATLAB

    % Example of optimization maximizing intensity in the target
    %
    % Copyright (C) 2019 Guilherme B Saturnino
    
    
    opt = opt_struct('TDCSoptimize');
    opt.leadfield_hdf = 'leadfield/ernie_leadfield_EEG10-10_UI_Jurak_2007.hdf5';
    opt.name = 'optimization/single_target';
    
    opt.max_total_current = 2e-3;
    opt.max_individual_current = 1e-3;
    opt.max_active_electrodes = 8;
    
    opt.target.positions = [-55.4, -20.7, 73.4];
    % Set the intensity to a large value (e.g, 100)
    opt.target.intensity = 100;
    
    % Run optimization
    run_simnibs(opt);
    

Using MNI Coordinates

The target positions are, as always in SimNIBS, given in world coordinates in subject space (see here for more information). However, we can use the mni2subject_coords function to transform coordinates from MNI space to subject space. When the transformed coordinates are outside the gray matter of the subject, they will be projected to the closest gray matter position.

  • Python

    import simnibs
    
    opt = simnibs.opt_struct.TDCSoptimize()
    opt.leadfield_hdf = 'leadfield/ernie_leadfield_EEG10-10_UI_Jurak_2007.hdf5'
    opt.name = 'optimization/MNI_target'
    
    opt.max_total_current = 2e-3
    opt.max_individual_current = 1e-3
    opt.max_active_electrodes = 8
    
    target = opt.add_target()
    # Transfrorm a set of coordinates from MNI space to subject space.
    # The second argument of the mni2subject_coords function
    # is the path to the "m2m_subID" folder.
    target.positions = simnibs.mni2subject_coords([-37, -21, 58], 'm2m_ernie')
    target.intensity = 0.2
    
    # Run optimization
    simnibs.run_simnibs(opt)
    

  • MATLAB

    opt = opt_struct('TDCSoptimize');
    opt.leadfield_hdf = 'leadfield/ernie_leadfield_EEG10-10_UI_Jurak_2007.hdf5';
    opt.name = 'optimization/MNI_target';
    
    opt.max_total_current = 2e-3;
    opt.max_individual_current = 1e-3;
    opt.max_active_electrodes = 8;
    
    % Transfrorm a set of coordinates from MNI space to subject space.
    % The second argument of the mni2subject_coords function
    % is the path to the "m2m_subID" folder.
    opt.target.positions = mni2subject_coords([-37, -21, 58], 'm2m_ernie');
    opt.target.intensity = 0.2;
    
    % Run optimization
    run_simnibs(opt);
    

Multiple targets

To optimize multiple distant targets simultaneously, just use multiple target structures.

  • Python

    ''' Example of an optimization with two targets
    '''
    
    import simnibs
    
    opt = simnibs.opt_struct.TDCSoptimize()
    opt.leadfield_hdf = 'leadfield/ernie_leadfield_EEG10-10_UI_Jurak_2007.hdf5'
    opt.name = 'optimization/two_targets'
    
    opt.max_total_current = 2e-3
    opt.max_individual_current = 1e-3
    opt.max_active_electrodes = 8
    
    # Target in the left motor cortex
    target_left = opt.add_target()
    target_left.positions = [-34.0, -21.4, 88.5]
    target_left.intensity = 0.2
    # Target in the right motor cortex
    target_right = opt.add_target()
    target_right.positions = [32.4, -25.5, 90.4]
    target_right.intensity = -0.2 # negative value revert the direction
    
    simnibs.run_simnibs(opt)
    

  • MATLAB

    % Example of an optimization with two targets
    %
    
    opt = opt_struct('TDCSoptimize');
    opt.leadfield_hdf = 'leadfield/ernie_leadfield_EEG10-10_UI_Jurak_2007.hdf5';
    opt.name = 'optimization/two_targets';
    
    opt.max_total_current = 4e-3;
    opt.max_individual_current = 2e-3;
    opt.max_active_electrodes = 8;
    
    % Target in the left motor cortex
    opt.target(1).positions = [-34.0 -21.4 88.5];
    opt.target(1).intensity = 0.2;
    % Target in the right motor cortex
    opt.target(2).positions  = [32.4 -25.5 90.4];
    opt.target(2).intensity = -0.2;
    
    
    run_simnibs(opt);
    

By using multiple targets, SimNIBS will try to hit each target with its intensity, whereas setting many positions in a single target, SimNIBS will try to hit the average intensity over the many positions.

Avoidance Regions

You can also add regions where the electric field should be more penalized than elsewhere. This is done using the avoid optional structure. In this examples, we will set the field to avoid the eyes.

  • Python

    ''' Example of an optimization punishing more the field in the eyes
    
        Copyright (C) 2019 Guilherme B Saturnino
    '''
    import simnibs
    
    opt = simnibs.opt_struct.TDCSoptimize()
    opt.leadfield_hdf = 'leadfield/ernie_leadfield_EEG10-10_UI_Jurak_2007.hdf5'
    opt.name = 'optimization/avoid'
    
    opt.max_total_current = 2e-3
    opt.max_individual_current = 1e-3
    opt.max_active_electrodes = 8
    
    target = opt.add_target()
    target.positions = simnibs.mni2subject_coords([-37, -21, 58], 'm2m_ernie')
    target.intensity = 0.2
    
    avoid = opt.add_avoid()
    avoid.tissues = 1006  # 1006 corresponds to the eye surface
    
    # Run optimization
    simnibs.run_simnibs(opt)
    

  • MATLAB

    % Example of an optimization punishing more the field in the eyes
    %
    % Copyright (C) 2019 Guilherme B Saturnino
    
    opt = opt_struct('TDCSoptimize');
    opt.leadfield_hdf = 'leadfield/ernie_leadfield_EEG10-10_UI_Jurak_2007.hdf5';
    opt.name = 'optimization/avoid';
    
    opt.max_total_current = 2e-3;
    opt.max_individual_current = 1e-3;
    opt.max_active_electrodes = 8;
    
    opt.target.positions = mni2subject_coords([-37, -21, 58], 'm2m_ernie');
    opt.target.intensity = 0.2;
    
    opt.avoid.tissues = 1006; % 1006 corresponds to the eye surface
    
    % Run optimization
    run_simnibs(opt)
    

Note

For more options and information on avoidance regions please see the referece for the TDCSavoid structure. You can visualize the position of the avoided region in the results by deselecting normE in gmsh, and selecting avoid_1.