EFMouse: 4x1 montage

EFMouse is a Matlab tool for electric field modelling in the mouse brain.

This notebook reproduces results for the 4x1 montage in:

Sanchez-Romero R., Akyuz, S., & Krekelberg, B. (2024). EFMouse: a Matlab toolbox to model electric fields in the mouse brain. bioRxiv. https://doi.org/10.1101/2024.07.25.605227

The notebook is also an introductory tutorial to simulate other montages.

(Developed by Ruben Sanchez-Romero and Bart Krekelberg,

Rutgers-Newark, Center for Molecular and Behavioral Neuroscience,

for support open an issue in https://github.com/klabhub/EFMouse/issues)

Table of Contents

Introduction Typical workflow Stage 0: Initialize Stage 1: Create electrodes and a craniotomy Stage 2: Find the electrode boundaries and add to the mesh Stage 3: export the mesh and model Stage 4: run getDP Tissue Based Analysis ROI Analysis Stage 5: Volumetric Analysis Atlas Based Analysis Tips and Tricks

Introduction

EFMouse combines several elements from previous work:

Anatomical data from the Digimouse project.
A refined finite element mesh of the Digimouse by Alekseichuk et al., (2019), Neuroimage.
Matlab code and inspiration from the ROAST toolbox.
The getDP finite element method (FEM) solver (v 3.2.0).
The Allen mouse brain atlas.

Key new functionality includes:

Positioning electrodes in arbitrary locations on the mouse body.
Modeling craniotomies in the mouse skull.
Analyzing the electric field in user-defined regions of interest
Analyzing the electric field in brain regions defined in the Allen mouse atlas.
A measure of relative focality of the electric field.
A measure of direction homogeneity of the electric field.

For details see Sanchez-Romero et al., (2024) and EFMouse.m Matlab code.

Typical workflow

The typical workflow to generate a model consists of several stages

stage 0 : Initialize.
stage 1: Define elecrodes and (optionally) a craniotomy.
stage 2: Impose boundary conditions.
stage 3: Export files for the FEM solver getDP.
stage 4: Run the FEM solver.
stage 5: Map the results to the Allen mouse atlas.

Stage 0: Initialize

We start by setting up the basic parameters of the EFMouse object:

o = EFMouse; % Create a default empty object of the class EFMouse
o.dir = '/Users/rubensanchez/desktop/EFMouse/4x1Montage'; % Results and the object (4x1.mat) will be saved here.
o.ID = '4x1';  % A name/tag for this simulation.
o.log  = true;  % Create a log file.

Note: The log file is created but doesn't update when running a simulation from a notebook (like here). It will properly update when running directly from the Matlab command line.

initialize() creates the folder to store the results. With overwrite=true, anything in this folder will be deleted.

o.initialize(overwrite=true);

The initialized EFMouse object contains the default mouse mesh (without stimulation electrodes). Let's visualize it.

plotMesh(o)

The mesh is pretty dense, so it will look like a black blob. The olfactory bulbs are visible on the top. In a regular figure window, you can use the Matlab figure tools to rotate or zoom. (You can also adapt the plotMesh() function in EFMouse.m to modifiy the figure directly.)

Stage 1: Create electrodes and a craniotomy

Now we define a 4x1 stimulation montage: 4 electrodes on the head (targeting visual areas), with a return electrode on the lumbar area.

o.eTag = ["Anterior" "Posterior" "Lateral" "Medial" "Lumbar"];
% currents need to sum to 0. (In milliamps (mA) units.)
o.eCurrent = [0.05,0.05,0.05,0.05,-0.2];
% center coordinate of the electrode position (in mesh space) 
o.eCenter = [-3.56,29.5,5.45;
             -3.43,24.8,6.02;
             -5.5,26.72,3.58; 
              -1.2,27.11,6.29;
              -1.1819,-9.655,10.129]';
% radius of the electrodes (in mesh space units) (1 mesh unit ~ 1 mm)
o.eRadius = [0.71,0.67,0.64,0.6,0.8]
o = EF Mouse Model (label: 4x1) in directory /Users/rubensanchez/desktop/EFMouse/4x1Montage (stage=0).

Add a craniotomy:

% center coordinate of the craniotomy position (in mesh space)
o.cCenter = [-3.4236,27.1067,5]';
% radius of the creaniotomy (in mesh space units) (1 mesh unit ~ 1 mm)
o.cRadius = 1.5;

To add these to the model, we run the pipeline to stage 1. With show=true, this also opens a figure for inspection of the mesh:

o.run(targetStage=1,show=true)

Compute stage 1 - computeMesh
EF Mouse Model (label: 4x1) in directory /Users/rubensanchez/desktop/EFMouse/4x1Montage (stage=0).
-Creating craniotomy-
 Craniotomy: touching tissue:
   gray: num elements = 6884
   csf: num elements = 6083
   bone: num elements = 8247
   skin: num elements = 1551
 Craniotomy only removes skin and bone.
-Creating 5 electrodes-
 Electrode: Anterior: touching tissue:
   skin: num elements = 39
   Creating electrode Anterior in max touched tissue: skin.
 Electrode: Posterior: touching tissue:
   skin: num elements = 46
   Creating electrode Posterior in max touched tissue: skin.
 Electrode: Lateral: touching tissue:
   skin: num elements = 125
   Creating electrode Lateral in max touched tissue: skin.
 Electrode: Medial: touching tissue:
   skin: num elements = 110
   Creating electrode Medial in max touched tissue: skin.
 Electrode: Lumbar: touching tissue:
   skin: num elements = 45
   Creating electrode Lumbar in max touched tissue: skin.

Stage 1 complete - 61.5275 seconds 

This mesh now shows the electrodes in red and blue (the blue return is out of the field of view; on the lumbar area), for the craniotomy, skin removal is shown in green and skull removal shown in yellow.

Stage 2: Find the electrode boundaries and add to the mesh

In this stage, EFMouse searches for the edges of the electrodes as this is where the current flows into the model. The boundary conditions of the Laplace equation are set here.

o.run(targetStage=2)
Compute stage 2 - computeBoundary
                 Electrode Area
                 ______________

    Anterior        0.78309    
    Posterior       0.67532    
    Lateral         0.59597    
    Medial          0.62975    
    Lumbar           1.0771    

Stage 2 complete - 85.6558 seconds 

Stage 3: export the mesh and model

Up until this stage, all changes were internal to the Matlab object (saved in 4x1.mat), now we export files that the FEM solver getDP can read. The .msh file contains information on the mesh (all the nodes, elements, labels and boundaries), and the .pro file tells getDP which partial differential equations (PDE) model it needs to solve this mesh.

o.run(targetStage=3)
Compute stage 3 - export data
----Starting saveMesh...29-Jan-2025 18:31:37
----saveMesh elapsed time: 34.8313 seconds

----Starting savePro...29-Jan-2025 18:32:12
----savePro elapsed time: 0.0565 seconds

Stage 3 complete - 53.6890 seconds 

For troubleshooting, have a look at the .pro file in the project folder; it contains all the model specifications, including conductivity (in siemens per meter S/m) for the different tissues/elements.

type(file(o,"PRO"))
/* 
 .pro file created by EFMouse on 29-Jan-2025 18:32:12 
 ID: 4x1 
 Dir: /Users/rubensanchez/desktop/EFMouse/4x1Montage 
 */

Group {

gray = Region[1];
csf = Region[2];
bone = Region[3];
skin = Region[4];
eye = Region[5];
craniotomy = Region[6];
skinremoved = Region[7];
Anterior = Region[8];
Posterior = Region[9];
Lateral = Region[10];
Medial = Region[11];
Lumbar = Region[12];
boundaryAnterior = Region[13];
boundaryPosterior = Region[14];
boundaryLateral = Region[15];
boundaryMedial = Region[16];
boundaryLumbar = Region[17];
DomainC = Region[{gray,csf,bone,skin,eye,craniotomy,skinremoved,Anterior,Posterior,Lateral,Medial,Lumbar}];
AllDomain = Region[{gray,csf,bone,skin,eye,craniotomy,skinremoved,Anterior,Posterior,Lateral,Medial,Lumbar,boundaryAnterior,boundaryPosterior,boundaryLateral,boundaryMedial,boundaryLumbar}];
}

Function {

sigma[gray] = 0.275;
sigma[csf] = 1.654;
sigma[bone] = 0.01;
sigma[skin] = 0.465;
sigma[eye] = 0.5;
sigma[craniotomy] = 1.654;
sigma[skinremoved] = 2.5e-14;
sigma[Anterior] = 5.9e+07;
sigma[Posterior] = 5.9e+07;
sigma[Lateral] = 5.9e+07;
sigma[Medial] = 5.9e+07;
sigma[Lumbar] = 5.9e+07;
du_dn1[] = 63.849704;
du_dn2[] = 74.039412;
du_dn3[] = 83.896792;
du_dn4[] = 79.396452;
du_dn5[] = -185.680706;
}

Jacobian {
  { Name Vol ;
    Case {
      { Region All ; Jacobian Vol ; }
    }
  }
  { Name Sur ;
    Case {
      { Region All ; Jacobian Sur ; }
    }
  }
}


Integration {
  { Name GradGrad ;
    Case { {Type Gauss ;
            Case { { GeoElement Triangle    ; NumberOfPoints  3 ; }
                   { GeoElement Quadrangle  ; NumberOfPoints  4 ; }
                   { GeoElement Tetrahedron ; NumberOfPoints  4 ; }
                   { GeoElement Hexahedron  ; NumberOfPoints  6 ; }
                   { GeoElement Prism       ; NumberOfPoints  9 ; } }
           }
         }
  }
}


FunctionSpace {
  { Name Hgrad_v_Ele; Type Form0;
    BasisFunction {
      // v = v  s   ,  for all nodes
      //      n  n
      { Name sn; NameOfCoef vn; Function BF_Node;
        Support AllDomain; Entity NodesOf[ All ]; }
    }
  }
}


Formulation {
  { Name Electrostatics_v; Type FemEquation;
    Quantity {
      { Name v; Type Local; NameOfSpace Hgrad_v_Ele; }
    }
    Equation {
      Galerkin { [ sigma[] * Dof{d v} , {d v} ]; In DomainC; 
      Jacobian Vol; Integration GradGrad; }
      Galerkin{ [ -du_dn1[], {v} ]; In boundaryAnterior ;                 Jacobian Sur; Integration GradGrad;}
      Galerkin{ [ -du_dn2[], {v} ]; In boundaryPosterior ;                 Jacobian Sur; Integration GradGrad;}
      Galerkin{ [ -du_dn3[], {v} ]; In boundaryLateral ;                 Jacobian Sur; Integration GradGrad;}
      Galerkin{ [ -du_dn4[], {v} ]; In boundaryMedial ;                 Jacobian Sur; Integration GradGrad;}
      Galerkin{ [ -du_dn5[], {v} ]; In boundaryLumbar ;                 Jacobian Sur; Integration GradGrad;}
    }
  }
}


Resolution {
  { Name EleSta_v;
    System {
      { Name Sys_Ele; NameOfFormulation Electrostatics_v; }
    }
    Operation { 
      Generate[Sys_Ele]; Solve[Sys_Ele]; SaveSolution[Sys_Ele];
    }
  }
}


PostProcessing {
  { Name EleSta_v; NameOfFormulation Electrostatics_v;
    Quantity {
      { Name v; 
        Value { 
          Local { [ {v} ]; In AllDomain; Jacobian Vol; } 
        }
      }
      { Name e; 
        Value { 
          Local { [ -{d v} ]; In AllDomain; Jacobian Vol; }
        }
      }
    }
  }
}
PostOperation {
{ Name Map; NameOfPostProcessing EleSta_v;
   Operation {
     Print [ v, OnElementsOf DomainC, File "/Users/rubensanchez/desktop/EFMouse/4x1Montage/4x1_v.pos", Format NodeTable ];
     Print [ e, OnElementsOf DomainC, Smoothing, File "/Users/rubensanchez/desktop/EFMouse/4x1Montage/4x1_e.pos", Format NodeTable ];
   }
}


}

This file can be opened in the getDP gui to run it manually, but stage 4 runs it for you.

Stage 4: run getDP

It will take getDP between 15 and 30 minutes to compute the solutions of the Laplace equation (on 2024 Mac or Windows hardware). The results are saved in _e.pos and _v.pos files (in getDP format). In the object the results are saved as .field and .voltage.

o.run(targetStage=4,show=false);
Compute stage 4 - computeField
Info    : Running '/Users/rubensanchez/Desktop/roast/EFmouse_class/EFMouse/lib/getdp-3.2.0/bin/getdpMac /Users/rubensanchez/desktop/EFMouse/4x1Montage/4x1.pro -solve EleSta_v -msh /Users/rubensanchez/desktop/EFMouse/4x1Montage/4x1.msh -pos Map' [GetDP 3.2.0, 1 node, max. 1 thread]
Info    : Started (Wed Jan 29 18:55:49 2025, Wall = 0.0974591s, CPU = 0.076878s, Mem = 5.34766Mb)
Info    : Initializing Gmsh
[1m[31mError   : Unknown number option 'General.NativeFileChooser'[0m
Info    : Loading problem definition '/Users/rubensanchez/desktop/EFMouse/4x1Montage/4x1.pro'
Info    : Selected Resolution 'EleSta_v'
Info    : Loading Geometric data '/Users/rubensanchez/desktop/EFMouse/4x1Montage/4x1.msh'
Info    : System 'Sys_Ele' : Real
[34mP r e - P r o c e s s i n g . . .[0m
Info    : Treatment Formulation 'Electrostatics_v'
  0%    : Pre-processing                                          
 10%    : Pre-processing                                          
 20%    : Pre-processing                                          
 30%    : Pre-processing                                          
 40%    : Pre-processing                                          
 50%    : Pre-processing                                          
 60%    : Pre-processing                                          
 70%    : Pre-processing                                          
 80%    : Pre-processing                                          
 90%    : Pre-processing                                          
                                          
Info    : System 1/1: 1027341 Dofs
Info    : (Wall = 47.7817s, CPU = 34.9588s, Mem = 683.68Mb)
[34mE n d   P r e - P r o c e s s i n g[0m
[34mP r o c e s s i n g . . .[0m
Info    : Generate[Sys_Ele]
  0%    : Processing (Generate)                                          
 10%    : Processing (Generate)                                          
 20%    : Processing (Generate)                                          
 30%    : Processing (Generate)                                          
 40%    : Processing (Generate)                                          
 50%    : Processing (Generate)                                          
 60%    : Processing (Generate)                                          
 70%    : Processing (Generate)                                          
 80%    : Processing (Generate)                                          
 90%    : Processing (Generate)                                          
                                          
Info    : Solve[Sys_Ele]
Info    : N: 1027341 - preonly lu mumps
Info    : SaveSolution[Sys_Ele]
Info    : (Wall = 1785.06s, CPU = 2042.32s, Mem = 5467.12Mb)
[34mE n d   P r o c e s s i n g[0m
[34mP o s t - P r o c e s s i n g . . .[0m
Info    : NameOfSystem not set in PostProcessing: selected 'Sys_Ele'
Info    : Selected PostProcessing 'EleSta_v'
Info    : Selected Mesh '/Users/rubensanchez/desktop/EFMouse/4x1Montage/4x1.msh'
Info    : PostOperation 1/2 
          > '/Users/rubensanchez/desktop/EFMouse/4x1Montage/4x1_v.pos'
  0%    : Post-processing (Compute)                                          
 10%    : Post-processing (Compute)                                          
 20%    : Post-processing (Compute)                                          
 30%    : Post-processing (Compute)                                          
 40%    : Post-processing (Compute)                                          
 50%    : Post-processing (Compute)                                          
 60%    : Post-processing (Compute)                                          
 70%    : Post-processing (Compute)                                          
 80%    : Post-processing (Compute)                                          
 90%    : Post-processing (Compute)                                          
                                          
Info    : PostOperation 2/2 
          > '/Users/rubensanchez/desktop/EFMouse/4x1Montage/4x1_e.pos'
  0%    : Post-processing (Generate)                                          
 10%    : Post-processing (Generate)                                          
 20%    : Post-processing (Generate)                                          
 30%    : Post-processing (Generate)                                          
 40%    : Post-processing (Generate)                                          
 50%    : Post-processing (Generate)                                          
 60%    : Post-processing (Generate)                                          
 70%    : Post-processing (Generate)                                          
 80%    : Post-processing (Generate)                                          
 90%    : Post-processing (Generate)                                          
                                          
  0%    : Post-processing (Compute)                                          
Info    : Smoothing (phase 1)
Info    : Smoothing (phase 2)
Info    : (Wall = 2072.26s, CPU = 2239.55s, Mem = 6005.12Mb)
[34mE n d   P o s t - P r o c e s s i n g[0m
Info    : Stopped (Wed Jan 29 19:30:24 2025, Wall = 2074.84s, CPU = 2239.98s, Mem = 6005.12Mb)
Stage 4 complete - 2158.7607 seconds 

Use the plotEf() function to visualize the electric field in the X direction, and then the electric field magnitude.

plotEf(o,type='eX',percentile=98, tissue='gray');

----Starting plotEf...29-Jan-2025 19:57:56

----plotEf elapsed time: 5.8017 seconds

plotEf(o,type='eMag',percentile=98,tissue='gray')

----Starting plotEf...29-Jan-2025 19:58:02

----plotEf elapsed time: 2.1372 seconds

Tissue Based Analysis

In general, we focus on the brain, but if necessary, we can compute electric field estimates for the rest of the mouse body for a full characterization of the anatomical effects of the electrical stimulation protocol. See Sanchez-Romero et al., (2024) for a full description of the electric field metrics.

analyzeTissue(o)
----Starting analyzeTissue...29-Jan-2025 19:58:39
Electric field summary statistics for gray tissue
              mean       median      std        min       max  
            ________    ________    ______    _______    ______

    eX        1.2703     0.99226    1.1544    -5.4246    14.526
    eY      -0.54963    -0.31191    1.5478    -12.784    12.106
    eZ       -1.6616     -1.2725    1.3786    -16.308    2.5609
    eMag      2.6406      2.1451    1.8255    0.21219    18.899

[Homogeneity ranges from 0 to 1]
Homogeneity = 0.8260, for ef_norm_mean: (0.507 -0.159 -0.632)

Electric field summary statistics for csf tissue
             mean       median      std        min       max  
            _______    ________    ______    _______    ______

    eX        1.794     0.89465    4.5906    -39.173      43.1
    eY      -1.8382    -0.61238    4.3158    -53.131    13.656
    eZ      -4.0961     -1.3658    8.0954    -95.879    35.001
    eMag     7.4538      4.1059     8.547    0.21219    102.53

[Homogeneity ranges from 0 to 1]
Homogeneity = 0.6418, for ef_norm_mean: (0.345 -0.151 -0.520)

Electric field summary statistics for bone tissue
             mean        median       std         min         max  
            _______    __________    ______    __________    ______

    eX      0.92165      0.020804    4.2643        -54.62    101.01
    eY      -1.8937       -1.2447    3.9041       -70.841    30.039
    eZ      -1.9948    -0.0036406    8.2483       -143.11    64.906
    eMag     5.1211        2.2468     9.146    4.0102e-06       151

[Homogeneity ranges from 0 to 1]
Homogeneity = 0.4253, for ef_norm_mean: (0.106 -0.409 -0.047)

Electric field summary statistics for skin tissue
              mean       median       std         min         max  
            ________    _________    ______    __________    ______

    eX       0.24094    0.0058755    2.2351        -105.3    114.15
    eY       -1.1117      -0.7269      2.28        -100.4    101.54
    eZ      -0.27302     0.015851    3.3009       -112.12    79.406
    eMag      1.9861       1.2332    4.3028    9.5456e-07    120.68

[Homogeneity ranges from 0 to 1]
Homogeneity = 0.5627, for ef_norm_mean: (0.060 -0.547 0.117)

Electric field summary statistics for eye tissue
              mean       median       std         min         max   
            ________    ________    ________    ________    ________

    eX       0.27523     0.29155    0.097919     0.10035     0.46299
    eY        0.1241    -0.22774     0.48241    -0.49078     0.95888
    eZ      -0.75129    -0.51244     0.36828     -1.5235    -0.31149
    eMag     0.95566     0.79114     0.34646     0.50607      1.7733

[Homogeneity ranges from 0 to 1]
Homogeneity = 0.8290, for ef_norm_mean: (0.355 -0.033 -0.748)

Electric field summary statistics for craniotomy tissue
             mean      median      std        min       max  
            _______    _______    ______    _______    ______

    eX       2.0635     1.2586    4.4655     -42.56    29.185
    eY      -2.8514    -2.2702    4.9046    -57.525    24.124
    eZ       -5.936    -3.1851    8.0263    -76.415    11.012
    eMag     8.2106     4.6355    9.4138     1.6988    87.469

[Homogeneity ranges from 0 to 1]
Homogeneity = 0.8691, for ef_norm_mean: (0.260 -0.464 -0.687)

Electric field summary statistics for skinremoved tissue
             mean      median      std        min       max  
            _______    _______    ______    _______    ______

    eX       3.4673     3.0961    8.0299     -42.56    29.185
    eY      -4.5407    -2.5214    9.4863    -52.951    24.124
    eZ      -14.601    -9.6814    14.492    -76.415    11.012
    eMag     19.066     13.873    15.704     1.5744    87.469

[Homogeneity ranges from 0 to 1]
Homogeneity = 0.8051, for ef_norm_mean: (0.312 -0.225 -0.707)

ROI Analysis

Now that we have the field estimates for each node in the mesh, we can do an ROI based analysis. For instance, a box ROI.

For the relative focality the reference area is the rest of the gray tissue.

roi.shape = 'box';

roi.dim = [[-4.0087 -2.56];[26.4022 27.64];[3 4]]';% Define the limits of the box ROI (x,y,z)(in mesh space units) (1 mesh unit ~ 1 mm)

summary = analyzeRoi(o,roi,plot=true,foc_percentile_max=99.9,foc_threshold=75); % Define values for relative focality metric.

----Starting analyzeRoi...29-Jan-2025 19:59:15
Electric field summary statistics for a box roi in gray
              mean       median       std        min       max  
            ________    ________    _______    _______    ______

    eX        3.0435      2.8444     1.1362    0.47771    8.5696
    eY      -0.47302    -0.35704    0.95461    -4.2281    1.6231
    eZ       -7.3414     -7.3198     1.2901    -10.475    -3.544
    eMag      8.0868      8.0123     1.3604     4.4533    12.903

[Relative focality ranges from 0 to 1]
Relative focality = 0.9895, with 161805 reference nodes (cutoff: eMag > 75.00% of the target area max (99.90th percentile))
[Homogeneity ranges from 0 to 1]
Homogeneity = 0.9855, for ef_norm_mean: (0.379 -0.061 -0.908)

% Rotate and zoom to view the (purple) box ROI underneath the craniotomy

xlim([-3.31 1.13])

ylim([26.3 36.1])

zlim([2.31 4.98])

view([62.61 27.62])

For comparison, define an homologous box ROI but in the right hemisphere

roi_control.shape = 'box';

roi_control.dim = [[2.56 4.0087];[26.4022 27.64];[3 4]]';% Define the limits of the box ROI (in mesh space units)

summary = analyzeRoi(o,roi_control,plot=true,foc_percentile_max=99.9,foc_threshold=75);

----Starting analyzeRoi...29-Jan-2025 19:59:24
Electric field summary statistics for a box roi in gray
              mean       median       std        min        max   
            ________    ________    _______    _______    ________

    eX        1.2782      1.2861    0.21287    0.70558      1.7305
    eY      -0.76965    -0.72403    0.14409    -1.1859    -0.53497
    eZ       -1.0586     -1.0639    0.11413    -1.3238    -0.73842
    eMag      1.8405      1.8344    0.19583     1.3633      2.2942

[Relative focality ranges from 0 to 1]
Relative focality = 0.3564, with 161865 reference nodes (cutoff: eMag > 75.00% of the target area max (99.90th percentile))
[Homogeneity ranges from 0 to 1]
Homogeneity = 0.9939, for ef_norm_mean: (0.692 -0.418 -0.579)

% Rotate and zoom to view the (purple) box ROI underneath the craniotomy

xlim([-3.31 1.13])

ylim([26.3 36.1])

zlim([2.31 4.98])

view([62.61 27.62])

Stage 5: Volumetric Analysis

The mesh coordinates are not particularly intuitive, and you may want to estimate electric fields in specific brain areas (as defined in an atlas). EFMouse does this in reference to the Allen Mouse Brain Atlas.

To use this, we first have to map the mesh-based results to the Allen Atlas. This is Stage 5 of the pipeline. This stage first exports the mesh-based results to a volume (using linear interpolation) and then uses the FLIRT tool in FSL to transform this volume to the coordinates of the Allen Atlas. Note that the alignment between the Digimouse mesh and the Allen Atlas is not perfect because they are based on different imaging modalities (and different mouse strains). (See Sanchez-Romero et al., (2014).)

This will fail if FSL is not installed.

o.run(targetStage=5,startStage=5)
Compute stage 5 - computeVoxelSpace
----Starting computeVoxelSpace...29-Jan-2025 20:09:14
----Exporting to NIFTI volumes 
---/Users/rubensanchez/desktop/EFMouse/4x1Montage/4x1_efm.nii.gz created
---/Users/rubensanchez/desktop/EFMouse/4x1Montage/4x1_efX.nii.gz created
---/Users/rubensanchez/desktop/EFMouse/4x1Montage/4x1_efY.nii.gz created
---/Users/rubensanchez/desktop/EFMouse/4x1Montage/4x1_efZ.nii.gz created
----computeVoxelSpace elapsed time: 40.9412 seconds

Atlas Based Analysis

Once stage 5 has completed we can query electric fields based on a region that is defined in the Allen Atlas. For instance: 'Visual areas' for the left hemisphere.

For the relative focality the reference area is the rest of the Allen atlas "Isocortex".

T = analyzeAtlas(o,"Visual areas",hemisphere="left",foc_threshold=75,foc_percentile_max=99.9,foc_reference='Isocortex');
----Starting analyzeAtlas...29-Jan-2025 20:21:05
Area: Visual areas (1.0% of brain) , hemisphere left 
              mean       median      std        min        max  
            ________    ________    ______    ________    ______

    eX        4.8488       4.937    2.5216      -2.194    12.861
    eY      -0.40598    -0.24019    2.0952     -10.013     7.306
    eZ       -6.3823      -6.751    2.5299     -13.498         0
    eMag       8.575      8.7693    2.8676    0.003294    16.273

[Relative focality ranges from 0 to 1]
Relative focality = 0.9985, with 92410 reference voxels (cutoff: eMag > 75.00% of the area max (99.90th percentile))
[Homogeneity ranges from 0 to 1]
Homogeneity = 0.9291, for ef_norm_mean: (0.558 -0.049 -0.741)
----analyzeAtlas elapsed time: 1.3862 seconds

For comparison, we can query results for the right hemisphere, which is contra-lateral from the targeted area.

T = analyzeAtlas(o,"Visual areas",hemisphere="right",foc_threshold=75,foc_percentile_max=99.9,foc_reference='Isocortex');
----Starting analyzeAtlas...29-Jan-2025 20:21:06
Area: Visual areas (1.0% of brain) , hemisphere right 
              mean       median       std         min         max  
            ________    ________    _______    _________    _______

    eX         1.706      1.8367    0.77312    0.0021152     3.0388
    eY      -0.76484    -0.82405    0.41637      -1.6543    0.14136
    eZ        -1.036     -1.1166    0.46465      -1.7651          0
    eMag       2.168      2.3652    0.92598    0.0025962     3.5464

[Relative focality ranges from 0 to 1]
Relative focality = 0.6212, with 92457 reference voxels (cutoff: eMag > 75.00% of the area max (99.90th percentile))
[Homogeneity ranges from 0 to 1]
Homogeneity = 0.9855, for ef_norm_mean: (0.782 -0.361 -0.479)
----analyzeAtlas elapsed time: 0.7801 seconds

Tips and Tricks

The run function can be called with any targetStage from 0 to 5. Stages that have already been completed will not be run again, and all stages that need to be completed to reach the targetStage will be completed in turn.
IMPORTANT: To continue working with a model from a previous Matlab session, you create an EFMouse object with the dir and ID input arguments to specify the folder with the (previously saved) results and the ID of the montage. For instance, to open the simulation from the tutorial and also know the stage:

o = EFMouse(dir='/Users/rubensanchez/desktop/EFMouse/4x1Montage',ID='4x1')