\title{FanPy Manual}
\author{
        Christian Schmidt \\
                cschmidt18057 at gmail.com
}
\date{\today}

\documentclass[12pt]{article}
\usepackage[latin1]{inputenc}
\usepackage[english]{babel}
\usepackage[T1]{fontenc}
\usepackage{hyperref}
\usepackage{graphicx}
\usepackage{siunitx}
\usepackage{listings}
\usepackage[a4paper]{geometry}
% Custom colors
\usepackage{color}
\definecolor{deepblue}{rgb}{0,0,0.5}
\definecolor{deepred}{rgb}{0.6,0,0}
\definecolor{deepgreen}{rgb}{0,0.5,0}
% Default fixed font does not support bold face
\DeclareFixedFont{\ttb}{T1}{txtt}{bx}{n}{10} % for bold
\DeclareFixedFont{\ttm}{T1}{txtt}{m}{n}{10}  % for normal

\begin{document}
\lstset{%
	frame=single,
	breaklines=true,
	keepspaces=true,
	otherkeywords={self},             % Add keywords here
	basicstyle=\ttm,
	keywordstyle=\ttb\color{deepblue},
	emphstyle=\ttb\color{deepred},    % Custom highlighting style
	stringstyle=\color{deepgreen}	
	}
\maketitle
\setcounter{tocdepth}{2} % Set the depth of the table of contents to show sections and subsections only

\newcommand{\fref}[1]{\mbox{(Fig. \ref{#1})}}

\tableofcontents % Print the table of contents


\section{Introduction}
FanPy (Field and Neuron Python) is a Python package for computing the field distribution during deep brain stimulation (DBS) and estimating the neural activation resulting from the applied electrical stimulus. Although the model pipeline is based on an application for deep brain stimulation, it is not limited to it and can be easily transfered to various bioelectrical applications. The Python package was used to estimate the neural activation extent in models of deep brain stimulation rendering different post-operative phases and target area properties, while reducing the computational expense by employing an approximation of the neural activation extent by field threshold values determined from the threshold-distance relationship of the applied axon models. For further details please refer to\newline
\textsf{Schmidt, C and van Rienen, U, Adaptive Estimation of the Neural Activation Extent during Deep Brain Stimulation, 2017, \href{https://arxiv.org/abs/1705.10478}{arXiv:1705.10478} [q-bio.NC]\newline}
The code is highly experimental and should be used with caution and a proper amount of healthy suspicion.\newline
FanPy is published under the GNU Lesser General Public License version 3.

\section{Tutorials}

\subsection{Generating a model of the subthalamic nucleus as target area}
\label{Modeling the Target Area}
\begin{figure}[!th]
    \centering
		\includegraphics[width=1.0\textwidth]{./pics/gm_mask_and_stn_model.PNG}
    \caption{Thresholding of the STN mask and the resulting surface model in 3D Slicer.}
    \label{fig:3dslicer}
\end{figure}
The subthalamic nucleus (STN) is a common target area during DBS. To incorporate the STN in a volume conductor model for DBS, the geometry information of the STN in an appropriate format, such as STEP, BREP, STL has to be available. If this geometry is not available, it has to be extracted from imaging data. Here we use the data available from the study\newline
\textsf{E. A. Accolla, J. Dukart, G. Helms, N. Weiskopf, F. Kherif, A. Lutti, R. Chowdhury, S. Hetzer, J. D. Haynes, A. A. K\"uhn, and B. Draganski, Brain tissue properties differentiate between motor and limbic basal ganglia circuits, Hum Brain Mapp, vol. 35, pp. 5083-5092, Oct 2014.}\newline
which is also available in the freely available software Lead DBS and provides a data set rendering the right and left STN and its sensorimotoric, limbic, and associative functional zones. The extent of the right and left STN are stored in \textsf{gm\_mask.nii}, which is processed with the open-source software 3D Slicer. First the right STN is extracted from the data using \textsf{Converts$\rightarrow$Crop Volume} and chosing a volume around it. Next the data is segmented by using \textsf{Segementation$\rightarrow$Editor} and applying the Threshold Effect with a range from $0.25$ to $1.00$ \fref{fig:3dslicer}. The resulting label map is used to build a surface model of the STN using \textsf{Surface Models$\rightarrow$Model Maker}. The resulting surface model is exported as an \textsf{.stl} (see \textsf{./files/STN.stl}).\newpage

\subsection{Generating a multi-compartment model geometry with Salome}
\label{Model Generation with Salome}
This tutorial shows the generation of DBS electrode lead model surrounded by multiple tissue compartments by using the open-source software \href{http://www.salome-platform.org/}{Salome}.

\paragraph{Salome Preferences}
The standard length unit in Salome is meters. Since the DBS electrode is in the scale of millimeters the precision in Salome has to be modified. An alternative would be to design the geomtry in Salome in meters and treating it as it would be in millimeters by scaling the mesh afterwards. In this tutorial the first approach is used, designing the electrode in millimeters in the Salome meter scale. For this, the precision under \textsf{File$\rightarrow$Preferences$\rightarrow$Geometry$\rightarrow$Input fields precision} has to be modified \fref{fig:salomepref}.
\begin{figure}[!th]
    \centering
		\includegraphics[width=1.0\textwidth]{./pics/salomepreferences.PNG}
    \caption{Preferences to be adjusted in SALOME in order to create the model in SI units.}
    \label{fig:salomepref}
\end{figure}

\paragraph{The DBS electrode lead model}
The DBS electrode lead model represents a \href{http://professional.medtronic.com/pt/neuro/dbs-md/prod/dbs-lead-model-3387}{Medtronic 3387 DBS electrode lead}, which comprises a cylindrical lead with a diameter of \SI{1.27}{mm} and four electrode contacts with a height of \SI{1.5}{mm} and a spacing of \SI{1.5}{mm}. The electrode will be positioned in a bounding box (representing the surrounding tissue) with an edge length of \SI{100}{mm} with the second electroce contact of the DBS electrode lead located in the center of the bounding box. We start be first creating a cylinder with height of \SI{1.5}{mm} and \SI{1.27}{mm} representing the second electrode contact \fref{fig:cylinder1}. The cylinder gets translated \SI{-0.75}{mm} along the z axis to be centered to the coordinate origin using \textsf{Transformation$\rightarrow$Translation}. Then we translate the cylinder again (which creates a copy of the cylinder) by \SI{-1.5}{mm} which represents the insulation between the second electrode contact and the first electrode contact. This is repeated to model the first electrode contact \fref{fig:cylinder12}.
\begin{figure}[!ht]
		\begin{minipage}[b]{.5\textwidth}
		\centering
		\includegraphics[width=1.0\textwidth]{./pics/cylinder1.PNG}
		\caption{Single cylinder}
		\label{fig:cylinder1}
	\end{minipage}
	\hfill
		\begin{minipage}[b]{.5\textwidth}
		\centering
		\includegraphics[width=1.0\textwidth]{./pics/cylinders12.PNG}
		\caption{Merged cylinders}
		\label{fig:cylinder12}
	\end{minipage}
\end{figure}\newline
The insulation below the first electrode contact has a rounded tip. We model this by creating a cylinder with a height of \SI{0.75}{mm} and half sphere with a diameter of \SI{1.27}{mm}. The half sphere is created by first creating a sphere and a cylinder with a height of \SI{0.635}{mm} (the radius of the sphere) and using \textsf{Operations$\rightarrow$Boolean$\rightarrow$Cut} to cut the cylinder out of the sphere, resulting in the half sphere. The cylinder and the half sphere are unified in one object using \textsf{Operations$\rightarrow$Boolean$\rightarrow$Fuse} and translated to the bottom of the first electrode contact \fref{fig:cylinder12tip}. Translating the second electrode contact cylinder multiple times, the third and forth electrode contact and their insulation are created \fref{fig:cylinder1234}.
\begin{figure}[!ht]
	\begin{minipage}[b]{.5\textwidth}
    \centering
		\includegraphics[width=1.0\textwidth]{./pics/cylinders12tip.PNG}
		\caption{Merged cylinders with electrode lead tip}
		\label{fig:cylinder12tip}
	\end{minipage}
	\hfill
	\begin{minipage}[b]{.5\textwidth}
    \centering
		\includegraphics[width=1.0\textwidth]{./pics/cylinders1234.PNG}
		\caption{Electrode contacts with electrode lead tip}
		\label{fig:cylinder1234}
	\end{minipage}
\end{figure}\newline
The lead of the electrode is a cylinder with a height of \SI{43.25}{mm} and positioned at the top of the fourth electrode contact. The seperated objects are then merged by using \textsf{New Entity$\rightarrow$Build$\rightarrow$Compound}. We still have duplicated faces at the intersection between two cylinders, which we can remove be using \textsf{Repair$\rightarrow$Glue Faces} with a tolerance of $1\cdot10^{-7}$ \fref{fig:electrode}. 
\begin{figure}[!th]
    \centering
		\includegraphics[width=1.0\textwidth]{./pics/electrode.PNG}
		\caption{DBS electrode geometry}
		\label{fig:electrode}
\end{figure}\newline
In the next step an encapsulation layer around the electrode with a thickness of \SI{0.2}{mm} is modeled using a cylinder with a diameter of \SI{1.67}{mm} and a height of \SI{54.5}{mm}. The half sphere at the encapsulation layer tip is created like that at the electrode tip, but using a diameter of \SI{0.835}{mm}. Finally, the cylinder and the half sphere are fused to form the encapsulation layer. Since the electrode geometry overlaps with the encapsulation layer object, they have to be seperated from each other in order to obtain an electrode geometry surrounded by an encapsulation layer geometry. For this the electrode geometry is cutted from the encapsulation layer object, creating a compound from the resulting object and the electrode geometry and finally glueing their faces using the functions applied in the previous steps \fref{fig:electrodewithencap}.
\begin{figure}[!th]
    \centering
		\includegraphics[width=1.0\textwidth]{./pics/electrodewithencap.PNG}
		\caption{DBS electrode with encapsulation layer geometry}
		\label{fig:electrodewithencap}
\end{figure}
The target area is represented by a model of the subthalamic nucleus, which modeling is described in the tutorial \ref{Modeling the Target Area}. The resulting \textsf{.stl} file is loaded into Salome and scaled using \textsf{Operation$\rightarrow$Transformation$\rightarrow$Scale Transform} by a factor of $0.001$ to transform the STN into the Salome meter coordinate system. The surface model of the STN is transformed into a solid using \textsf{New Entity$\rightarrow$Build$\rightarrow$Solid} and translated with the negative coordinates of the center of the motor region to position the STN at the second active electrode contact (which is in the origin of the coordinate system) \fref{fig:stngeometry}. To build the distinct compartments of the electrode, the encapsulation layer, the STN, and the bounding box, the geometries have to be separated from each other. First, the encapsulation layer object is cutted from the STN object resulting in the STN geometry. Then the electrode model with the encapsulation layer is partitioned using the STN geometry using \textsf{Opterations$\rightarrow$Partition}. This is required to generated common edges and faces at the intersection between STN and encapsulation layer. In the next step, the electrode geometry including the encapsulation layer and the STN geometry are merged into a compound and duplicate faces are removed by glueing their common faces \fref{fig:stngeometrydone}.
\begin{figure}[!ht]
	\begin{minipage}[b]{.5\textwidth}
    \centering
		\includegraphics[width=1.0\textwidth]{./pics/stn_imported.PNG}
		\caption{Imported STN geometry, which is not water-tight connected to the electrode geometry}
		\label{fig:stngeometry}
	\end{minipage}
	\hfill
	\begin{minipage}[b]{.5\textwidth}
    \centering
		\includegraphics[width=1.0\textwidth]{./pics/stn_geometry.PNG}
		\caption{STN geometry is water-tight connected and glued to the electrode with encapsulation layer geometry.}
		\label{fig:stngeometrydone}
	\end{minipage}
\end{figure}\newline
In the last step the bounding box, representing the surrounding tissue is modeled by creating a Box (cube) with an edge length of \SI{100}{mm} and translating its center to the origin of the coordinate system. The electrode geometry including the encapsulation layer and the STN is cutted from the bounding box resulting in the tissue box geometry, which is merged to a compound together with the electrode-encapsulation-STN geometry. Finally their common faces are glued, which results in the final model geometry of the DBS electrode surrounded by an encapsulation layer, the STN as target area and the whole model geometry surrounded by a tissue bounding box. The model geometry is exported in BREP format (see \textsf{./files/dbs\_electrode\_with\_stn\_motor.brep}). The whole Salome model pipeline is saved in hdf5 format (see \textsf{./files/dbs\_electrode\_stn\_motor.hdf}) \fref{fig:modelgeometry}.
\begin{figure}[!th]
    \centering
		\includegraphics[width=1.0\textwidth]{./pics/full_dbs_electrode_model.PNG}
		\caption{Full model geometry of the DBS electrode, encapsulation layer, STN, and the surrounding bounding box.}
		\label{fig:modelgeometry}
\end{figure}\clearpage
\newpage

\subsection{Meshing the model geometry}
\label{Mesh Generation}
For the meshing of the geometry the freely available software \href{http://gmsh.info}{Gmsh} is used. The meshing of the geometry requires several steps, from identifying subdomains (model compartments) and special boundaries (compartment surfaces) on which later the material properties and boundary conditions will be applied, to manually refining mesh sizes in specific subdomains or on specific boundaries. For this purpose a \textsf{.geo} file is created in which the model geometry created in \ref{Model Generation with Salome} is imported as a \textsf{.brep} file:

\begin{lstlisting}
Merge "dbs_electrode_with_stn_motor.brep";
\end{lstlisting}
A way to identify the specific subdomains (volumes) of the model compartments and the special boundaries is the visibility tool under \textsf{Tools->Visibility}. In the \textsf{Tree browser} the model parts matching a certain Volume, Surface, or Point are shown when selected \fref{fig:gmsh}. For an easier use in the later field models, the information on the subdomains and boundaries can be written as comments in the \textsf{.geo} file.
\begin{figure}[!th]
    \centering
		\includegraphics[width=.9\textwidth]{./pics/gmsh_visibility.PNG}
		\caption{Setting certain model structures visible to identify the corresponding surfaces and points in Gmsh.}
		\label{fig:gmsh}
\end{figure}\newpage
\begin{lstlisting}
// STN: Volume 1
// Contact 4: Volume 9
// Contact 3: Volume 7
// Contact 2: Volume 2
// Contact 1: Volume 4
// Insulation Lead : Volume 10
// Insulation 4-3 : Volume 8
// Insulation 3-2 : Volume 6
// Insulation 2-1 : Volume 3
// Insulation 1-t : Volume 5
// Encapsulation : Volume 11
// Bounding Box : Volume 12
\end{lstlisting}
For the assignment of the material properties, it is not necessary to distinguish between model compartments with the same material properties. For example, the insulation between electrode contact 1 and 2 and electrode contact 2 and 3 does have the same material properties and can be merged to a so called \textsf{Physical Group}, which is a function provided by Gmsh to unify different subdomains to one label. \textbf{It is mandatory to label all sub-volumes of the model geometry with Physical Volumes, since only those get meshed by Gmsh.}
\begin{lstlisting}
// Boundaries
// Contact 1
Physical Surface(1) = {1386};
// Contact 2
Physical Surface(2) = {1381};
// Contact 3
Physical Surface(3) = {1392};
// Contact 4
Physical Surface(4) = {1396};
// Ground
Physical Surface(10) = {1404:1409};
\end{lstlisting}\newpage
\begin{lstlisting}
// Subdomains (Caution! Have to be complete, because only physical groups get meshed and exported by gmsh!)
// Contacts
Physical Volume(100) = {2, 4, 7, 9};
// Insulation
Physical Volume(101) = {3, 5, 6, 8, 10};
// Encapsulation
Physical Volume(102) = {11};
// Tissue
Physical Volume(103) = {12};
// STN
Physical Volume(104) = {1};
\end{lstlisting}
For the manual mesh refinement, the \textsf{Characteristic Length} of model points can be set in Gmsh, which determines the maximum element size for this point and neighbouring entities of higher dimension (lines, surfaces, volumes). For example, the surface of the electrode contacts, the encapsulation layer, and the target area should be meshed with a finer resolution, since the applied stimulation signal causes the largest changes in the field solution in this region. The surrounding tissue bounding box can be meshed with a coarser resolution. The \textsf{Visibility Tool} can be used again to identify these points. For this a certain surface or volume can be selected in the \textsf{Tree browser} (Model entities) which highlights the corresponding points as well \fref{fig:gmshselectpoints}.
\begin{figure}[!th]
    \centering
		\includegraphics[width=1.0\textwidth]{./pics/gmsh_select_points.png}
		\caption{Selecting surfaces in Gmsh to obtain the corresponding points for applying mesh refinement via \textsf{Characteristic Length}.}
		\label{fig:gmshselectpoints}
\end{figure}\newline
\begin{lstlisting}
// Characteristic Length for Bounding Box
Characteristic Length {744:754} = 5.0e-3;

// Characteristic Length for STN
Characteristic Length {1:732} = 0.2e-3;

// Characteristic Length for Encapsulation
Characteristic Length {28:29, 36:38, 40:87, 733:746} = 0.1e-3;

// Characteristic Length for Insulation
Characteristic Length {733:743} = 0.1e-3;

// Characteristic Length for Contacts
Characteristic Length {733:736, 739:742} = 0.1e-3;
\end{lstlisting}
In addition to setting the \textsf{Characteristic Length} at the desired points, it is also possible to define a region in which the mesh should be refined to a certain value.
\begin{lstlisting}
Field[1] = Box;
Field[1].VIn = 0.5e-3;
Field[1].VOut = 5e-3;
Field[1].XMax = 0.01;
Field[1].XMin = -0.01;
Field[1].YMax = 0.01;
Field[1].YMin = -0.01;
Field[1].ZMax = 0.01;
Field[1].ZMin = -0.01;
Background Field = 1;
\end{lstlisting}
Currently FEniCS can only import mesh elements of order 1. To restrict the meshing in Gmsh to this order and to further optimize the mesh quality the following parameters are set.
\begin{lstlisting}
Mesh.ElementOrder = 1;
Mesh.Optimize = 1;
\end{lstlisting}
All the information is stored in a \textsf{.geo} file (see \textsf{./files/dbs\_electrode\_with\_stn\_motor.geo}), which can be meshed with Gmsh via command line \fref{fig:mesh}.
\begin{figure}[!th]
    \centering
		\includegraphics[width=1.0\textwidth]{./pics/mesh.png}
		\caption{Resulting mesh of the DBS electrode model illustrating the manual mesh refinement at the electrode contacts, encapsulation layer, and the STN.}
		\label{fig:mesh}
\end{figure}

\subsection{Solving the field equation with FEnICS}
\label{Solving the Field Equation with FEniCS}
The Field Model Field Model in \textsf{FanPy} is currently able to compute the solution to the field equation for stationary current field problems (SF) and electro-quasistatic field problems (EQS). The stationary current field problem is based on solving the Laplace equation for heterogeneous, purely resistive material properties, i.e. defined by an electrical conductivity. The electro-quasistatic field problem is based on solving the frequency-domain Laplace equation for heterogeneous, complex-valued material properties, i.e. defined by an electrical conductivity and relative permittivity at a specific angular frequency. The boundary conditions comprise Dirichlet boundary conditions, i.e. used for setting the potential at the surface of the active electrode contact and ground at the exterior boundary surfaces of the bounding box, and Neumann boundary conditions, i.e. setting insulation on the surfaces of the insulating parts of the electrode lead. For further details refer to the paper referenced in the readme of the FanPy repository.
The mesh of the model geometry on which the field equation with boundary conditions is solved is generated from the created \textsf{.geo} file using Gmsh (see \ref{Mesh Generation}) and converting it into a format readable by FEniCS using the command line tool \textsf{dolfin-convert} included in the FEniCS installation.

\paragraph{Setting up the Model Parameters using a XML Configuration File}
To use a model geometry in the \textsf{.geo} format created with Salome and Gmsh with the Field Model in \textsf{FanPy}, a \textsf{.xml} configuration file is required, which serves, on the one hand, to determine the default model parameters required to solve the field equation for the corresponding model and, on the other hand, to validate the physical volumes and surfaces defined in the \textsf{.geo} file and the model parameters (material properties, boundary conditions) set in the \textsf{.xml} configuration file. For the DBS model created in \ref{Model Generation with Salome} and \ref{Mesh Generation}, an example \textsf{.xml} configuration file is given by:

\begin{lstlisting}[language=xml]
<fmodel>
    <geofile>dbs_electrode_with_stn_motor</geofile>
    <geofile_debug>dbs_electrode_with_stn_motor_debug</geofile_debug>
    <subdomains>
        <subdomain name="electrode_contacts" material="platinum-iridium" physid="100" group="electrode"/>
        <subdomain name="electrode_lead" material="insulation" physid="101" group="electrode"/>
        <subdomain name="encapsulation" material="encapsulation" physid="102"/>
        <subdomain name="brain" material ="brain" physid="103"/>
        <subdomain name="stn" material="stn" physid="104"/>
    </subdomains>
    <boundaries>
        <boundary name="contact1" physid="1"/>
        <boundary name="contact2" physid="2"/>    
        <boundary name="contact3" physid="3"/>
        <boundary name="contact4" physid="4"/>
        <boundary name="ground" physid="10"/>
    </boundaries>
    <materials>
        <material name="platinum-iridium" conductivity="1e7"/>
        <material name="insulation" conductivity="1e-7"/>
        <material name="encapsulation" conductivity="0.15"/>
        <material name="brain" conductivity="0.30"/>
        <material name="stn" conductivity="0.20"/>
    </materials>
    <boundary_conditions>
        <boundary_condition name="active_electrode" type="dirichlet" boundary="contact2" value="1.0"/>
        <boundary_condition name="ground" type="dirichlet" boundary="ground" value="0.0"/>
    </boundary_conditions>       
</fmodel>
\end{lstlisting}

\begin{itemize}
	\item \textsf{fmodel}: root node for the FEniCS model.
	\item \textsf{geofile}: name of the \textsf{.geo} file without file ending
	\item \textsf{geofile\_debug}: name of the \textsf{.geo} file for debug mode (e.g. coarser mesh version)
	\item \textsf{subdomains}: root note for the subdomain parameters (sub-volumes)
		\begin{itemize}
			\item \textsf{subdomain}: node for subdomain parameters
			\begin{itemize}
				\item \textsf{name}: attribute for the name of the subdomain
				\item \textsf{material}: attribute for the name of the material defined at `materials\} 
				\item \textsf{physid}: attribute for the id of the corresponding physical group as defined in the \textsf{.geo} file
				\item \textsf{group}: (optional) attribute for merging subdomains to a group
			\end{itemize}
		\end{itemize}
	\item \textsf{boundaries}: root note for the boundary parameters (surfaces)
		\begin{itemize}
			\item \textsf{boundary}: node for boundary parameters
			\begin{itemize}
				\item \textsf{name}: attribute for the name of the boundary
				\item \textsf{physid}: attribute for the id of the corresponding physical group as defined in the \textsf{.geo} file
			\end{itemize}
		\end{itemize}
	\item \textsf{materials}: root note for the materials and their default values
		\begin{itemize}
		\item \textsf{material}: node for the material parameters
			\begin{itemize}
				\item \textsf{name}: attribute for the name of the material
				\item \textsf{conductivity}: attribute for the default electrical conductivity value in S/m
				\item \textsf{permittivity}: (optional) attribute for the default relative permittivity value (\textbf{required for EQS field problems})
			\end{itemize}
		\end{itemize}
	\item \textsf{boundary\_conditions}: root note for the boundary conditions and their default value
		\begin{itemize}
		\item \textsf{boundary\_condition}: note for the boundary condition parameters
			\begin{itemize}
				\item \textsf{name}: attribute for the name of the boundary condition
				\item \textsf{type}: attribute for the type of the boundary condition (currently only dirichlet supported in configuration file)
				\item \textsf{boundary}: attribute for the name of the boundary defined at \textsf{boundaries}
				\item \textsf{value}: attribute for the default electric potential value in V
			\end{itemize}
		\end{itemize}
	\item \textsf{frequency}: (optional) frequency value for which field solution should be computed in Hz (\textbf{required for EQS field problems})
\end{itemize}

\paragraph{Implementing the Model in Python using the FanPy FieldModel}
The model can be implemented in the Python package by inheriting the Field Model with setting the path to the model \textsf{.xml} file under \textsf{model\_path} and the file name of the \textsf{.xml} file under \textsf{MODELNAME}.
\begin{lstlisting}[language=python]
class VolumeConductorDBSElectrodeWithSTNInMotorRegion(AbstractFenicsFieldModel):
       
    def __init__(self):
        AbstractDBSFieldModel.__init__(self)
        self.MODELNAME = 'dbs_electrode_with_stn_motor'
        self.model_path +='dbs_electrode_with_stn/'
\end{lstlisting}
To set all required variables and parameters to run the model, the model has to be initialized with
\begin{lstlisting}[language=python]
fieldmodel = VolumeConductorDBSElectrodeWithSTNInMotorRegion()
fieldmodel.init()
\end{lstlisting}
If the model should be run with the default model parameters for the material properties, an empty dictionary can be given as argument.
\begin{lstlisting}[language=python]
fieldmodel.run(dict())
\end{lstlisting}
If the material properties should be varied, the dictionary has to define the material properties according their names defined in the \textsf{.xml} model file. For varying the material properties and rerunning the model, the model does not have to be reinitialized.
\begin{lstlisting}[language=python]
fieldparams = {
               'dielectric_properties': {
                  'brain': {'conductivity': 0.3},
                  'encapsulation': {'conductivity': 0.15},
                  'stn': {'conductivity': 0.5}
                  },
              }
fieldmodel.run(fieldparams)
\end{lstlisting}
\newpage
\paragraph{Evaluating the Field Solution and Derived Quantities}
The computed field solutions can be evaluated at a set of coordinates in the computational domain. The set of coordinates should be in the format of a \textsf{(n,dim)} array, with \textsf{dim=3} for a 3D field problem. Currently supported field quantities are:

\begin{itemize}
	\item electric potential: \textsf{fieldmodel.get\_potential()}
	\item electric field norm: \textsf{fieldmodel.get\_electric\_field\_norm()}
	\item current density norm: \textsf{fieldmodel.get\_current\_density\_norm()}
\end{itemize}
The field quantity (here the electric potential) is then evaluated with the following line
\begin{lstlisting}[language=python]
potential = fieldmodel.evaluate_quantity(fieldmodel.get_potential(), coordinates) 
\end{lstlisting}
with \textsf{coordinates} being a \textsf{numpy (n,dim) array}.
In addition to field equantities, integral quantities can be computed with:
\begin{itemize}
	\item Total power: \textsf{fieldmodel.get\_total\_power()}
	\item Impedance: \textsf{fieldmodel.get\_impedance()}
\end{itemize}

\paragraph{Storing and Loading the Model State}
The solution of a fieldmodel can be stored using
\begin{lstlisting}[language=python]
fieldmodel.store_model_state(save_path, save_name)
\end{lstlisting}
with \textsf{save\_path} defining the path and \textsf{save\_name} defining the name without file-ending. A \textsf{.hdf5} file containing meta information on the model, such as model parameters and order of the basis functions, and \textsf{.hdf5} file containing the model data, such as mesh and solution, are created.\newline
The solution of a field model can be loaded using
\begin{lstlisting}[language=python]
fieldmodel.load_model_state(save_path, save_name)
\end{lstlisting}

\subsection{Analytical Validation of the field solution}
The example problem consists of determining the electric potential distribution and electric field norm distribution for a layered sphere located in an homogeneous electric field \fref{fig:exampleproblem}. Details on solving the stationary current field problem and electro-quasistatic field problem for this geometry can be found in the paper mentioned in the package readme.
\begin{figure}[!th]
    \centering
		\includegraphics[width=.5\textwidth]{./pics/analytical_problem_geometry.png}
		\caption{Example geometry of a layered sphere inside a homogeneous electric field.}
		\label{fig:exampleproblem}
\end{figure}
The geometry parameters for the computational model are set to the following:
\begin{table}[!t]
	\centering
  \renewcommand{\arraystretch}{1.3}
	\caption{\label{tab:exampleparams} Geometric parameters and electric properties.\newline*Only required for electro-quasistatic field problem.}
	\begin{tabular}{llc}
			\hline
			\bfseries Compartment & \bfseries Property 	& \bfseries Value \\
			\hline
			Inner Sphere 			& Radius									& \SI{5}{mm} \\
												& Conductivity						& \SI{2.0}{S/m}\\
												& Relative Permittivity* 	& $120$\\
			Outer Sphere 			& Radius 									& \SI{15}{mm}\\
												& Conductivity						& \SI{0.1}{S/m}\\
												& Relative Permittivity* 	& $2\cdot10^6$\\
			Bounding Box 			& Edge length							& \SI{100}{mm} \\
												& Conductivity 						& \SI{1.0}{S/m}\\
												& Relative Permittivity* 	& $80$\\
			\hline
												& Frequency*							& \SI{35}{kHz} \\
			\hline
	\end{tabular}	
\end{table}
A potential of 0.5 V and - 0.5 V was applied to the left and right surface of the bounding box, respectively, to generate a homogeoneous field of 10 V/m in the bounding box.
For solving the stationary current field problem the conductivities in Table \ref{tab:exampleparams} are used (for solving the electro-quasistatic field problem the additional model parameters for the relative permittivity and the corresponding frequency from Table \ref{tab:exampleparams} are applied). The example problem can be carried out with FanPy by using the \textsf{AnalyticFieldModel} class. An example study can be carried out in FanPyStudies using
\begin{lstlisting}[language=python]
AnalyticValidationModel().run()
AnalyticValidationModel().postprocessing()
\end{lstlisting}
The model files to be used with the can be used with \textsf{AbstractFenicsFieldModel} class can be found in the local files:
\begin{itemize}
	\item Salome Model Geometry: \textsf{./files/analytic\_example\_model.hdf}
	\item Exported Model Geometry: \textsf{./files/analytic\_example\_model.brep}
	\item Gmsh Geo File:	\textsf{./files/analytic\_example\_model.geo}
	\item XML Model File: \textsf{./files/analytic\_example\_model.xml}
\end{itemize}
 
\paragraph{Comparison of the Analytic and Numeric Field Solution}
The following image shows a comparison between the analytic and numeric field solution for the example model with the model parameters mentioned above \fref{fig:examplesolution}.
\begin{figure}[!th]
    \centering
		\includegraphics[width=1.0\textwidth]{./pics/result_analytic_example.png}
		\caption{Comparison of the numeric with the analytic field solution.}
		\label{fig:examplesolution}
\end{figure}\newpage

\subsection{Computing the neural activation extent during DBS}
During deep brain stimulation (DBS), the applied stimulation signal on the electrode contacts of the DBS electrode results in a stimulation of the neuronal tissue in the proximity of the stimulation site. To estimate the extent of neural activation in the target area, the computational approach of determining the activation of perpendicularly oriented axons positioned around the electrode is commonly applied (Volume of Tissue Activated: VTA). For further details on the estimation of the neural activation extent refer to the paper mentioned in the repository description.

\paragraph{DBS Signal}
The extent of neural activation depends on the applied stimulation signal. In DBS for movement disorders, such as Parkinson's disease, commonly a cathodic (negative stimulation voltage or current with respect to the equilibrium potential), rectangular pulse with a frequency of \SI{130}{Hz} and pulse duration of \SI{60}{\micro s} is applied. Common stimulation protocols comprise voltage-controlled and current-controlled stimulation, depending on if the stimulation amplitude during the pulse duration is a constant voltage or a constant current, respectively.

The time-dependent field solution of the field model (see \ref{Solving the Field Equation with FEniCS}) resulting from the applied stimulation signal is obtained by multiplying the field solution (e.g., the electric potential) at a coordinate with the signal. While this approach is applicable for stationary field problems and electro-quasistatic field problems with constant frequency, a dispersive field problem would require the convolution of the field solution with the signal by using the Fourier Finite Element Method. The field solution for a monopolar (only one active electrode contact and one ground) current-controlled stimulation signal is obtained in a similar way with the field solution mutiplied by the impedance of the model, which is equivalent to scale the unit potential applied to the active electrode contact by the impedance in order to obtain a unit current.\newline
In FanPy a DBS signal can be created by

\begin{lstlisting}[language=python]
signal_params = {
	DBSsignal.FREQUENCY: 130,
	DBSsignal.AMPLITUDE: -1e-3, # 1 mA
	DBSsignal.DT: 5e-6,
	DBSsignal.NUMBER_OF_PULSES: 10,
	DBSsignal.PULSE_DURATION: 6e-5}
dbs_signal = DBSsignal(signal_params) 
\end{lstlisting}
which represents a train of $10$ rectangular current-controlled DBS signals with a frequency of \SI{130}{Hz}, a stimulation amplitude of \SI{1}{mA}, a pulse duration of \SI{60}{\micro s}, and a time step of \SI{5}{\micro s}.

\paragraph{Axon Model}
The axon model used for determining the activation with respect to an applied time-dependent electric potential is based on the model of a mammalian nerve fiber published in \href{https://www.ncbi.nlm.nih.gov/pubmed/11826063}{Cameron McIntyre, et al, 2002}. The model is implemented in \href{https://www.neuron.yale.edu/neuron/}{NEURON}, which is an open-source software tool for modeling and simulating neurons, neuronal networks, their dynamics, and response to external stimuli.

The model comprises a multi-compartment double cable model of an axonal fiber with 21 nodes of Ranvier and 10 internodal segments resulting in a total of 221 compartments. The time-dependent electric potential along the axon is applied as an extracellular potential to each of the compartments in order to compute the inner potential of the axon model for detecting whether an action potential is elicited.

In FanPy, the computation of the activation of the axon model is carried out with the \textsf{Threshold\_Volume} class. For a set of axons added to the Threshold Volume, the minimum stimulation amplitude required to elicit an action potential for each axon model is computed. In the following example the threshold-distance relationship for $8$ axons with a diameter of \SI{5.7}{\micro m} starting in a distance of \SI{1}{mm} the active electrode contact center and a spatial stepping of \SI{0.5}{mm} should be computed. For computing the applied extracellular potential see \ref{Solving the Field Equation with FEniCS}.
\newpage
\begin{lstlisting}[language=python]
import numpy as np
import matplotlib.pyplot as plt

axons = Threshold_volume()
start_distance = 1e-3
spatial_step = 0.5e-3
ncpus = 1	# number of cpus, that should be used for computation
naxons = 8
for i in xrange(naxons):
	location = Threshold_Location(0,0,i)
	point = Threshold_Point(location)
	transform = point.get_transform()
	coordinate = np.array([0, start_distance+location.get_normal()*spatial_step, 0])	# axons are positioned on y-axis
	axon = Axon({'diameter': 5.7})	# 5.7 um axon model
	nodes = transform.transform_coords(axon.get_nodes())	# transform axon nodes to the defined coordinate
	point.set_coordinate(nodes[0]/2+nodes[-1]/2)
	point.set_axon(axon)
	point.set_dbs_signal(dbs_signal)
	potential = fieldmodel.evaluate_quantity(fieldmodel.get_potential(), nodes)
	point.set_potential(potential)
	axons.add_point(point)

axons.compute_optimal_activation(ncpus, 1.0)
ct_matrix = axons.get_coord_threshold_matrix()
unit_amplitude = dbs_signal.get_params()[DBSsignal.AMPLITUDE]
plt.plot(ct_matrix[:,1],ct_matrix[:,3]*unit_amplitude)	# plot y-axis of coordinates versus the stimulation amplitude
plt.xlabel('distance [m]')
plt.ylabel('Stimulation amplitude [mA]')
if unit_amplitude < 0:
	plt.gca().invert_yaxis()
plt.show()
\end{lstlisting}

\paragraph{Adaptive Estimation of the Neural Activation Extent}
\begin{figure}[!th]
    \centering
		\includegraphics[width=0.8\textwidth]{./pics/visualize_adaptive_algorithm.png}
		\caption{Illustration of the adaptive algorithm to compute the neural activation during DBS. Starting at seed points (e.g. placed in front of the active electrode contact) the algorithm determines whether an axon at this location would be activated by the applied stimulus. If yes, it places new seed points above, below, and right to the current seed point. This procedure continues until a closed zero-activation hull (seed points with no activation) is found. After removal of unneeded seed points to determine threshold isolines in the desired stimulation amplitude range, the minimum required stimulation amplitude to elicit an action potential in each axon model is determined via a bisectioning approach. From this threshold map, the VTA can be determined by deriving the ioslines (isosurfaces) for each stimulation amplitude within the given stimulation amplitude range.}
		\label{fig:adaptivealgorithm}
\end{figure}
Based on the model parameters, such as the dielectric properties of the tissue, the neural activation extent for the same DBS signal may vary, which requires pre-knowledge on how many axons have to be distributed in the target area in order to estimate the neural activation extent for a given stimulus. To overcome this pre-knowledge, an adaptive approach is implemented in FanPy, which uses seed points in the target area to estimate the region of activation for a given stimulation amplitude range \fref{fig:adaptivealgorithm}. First the activation region for the minimum stimulation amplitude is determined. Starting from the determined axon locations, the activation region fo the maximum stimulation amplitude is determined. In order to further save computation time, all points, which are not necessary to determine the neural activation extent within the stimulation amplitude range, are removed. For the remaining axon models the minimum stimulation amplitude required to elicit an action potential is computed, resulting in a threshold map. From this treshold map the neural activation extent is determined by computing isolines for a given stimulation amplitude within the stimulation amplitude range. The neural activation extent can be determined by assuming rotational symmetry, which requires only the estimation of the neural activation extent in one plane, or assuming rotational asymmetry, which computes the neural activation extent in `n` planes around the stimulation electrode. 
An example code in FanPy to estimate the neural activation extent for a given stimulation signal, stimulation amplitude range, and field solution with assuming rotational symmetry is given by

\begin{lstlisting}[language=python]
ncpus = 1	# number of cpus to be used for computation
neuronparams = {
		'g_step': 0.5e-3,			# spatial stepping
		'rot_degree_step': 10,		# rotational degree stepping
		'dbs_signal': dbs_signal,	# dbs signal with unit amplitude
		'fiberD': 5.7,				# axon fiber diameter
		'stim_min': 2.0,			# minimum stimulation threshold
		'stim_max': 4.0,			# maximum stimulation threshold
		'symmetric': True			# assume rotational symmetry
			}
vta = vtautil.compute_vta(ncpus, neuronparams)
\end{lstlisting}

Often it is of interest which stimulation amplitude range is required to obtain a neural activation extent within a specified distance. For example, in DBS of the subthalamic nucleus (STN), the stimulation amplitude, which activates a certain region in the (STN), but not exceeds a certain distance to reduce side effects, is required.
One approach to determine the required stimulation amplitude range is to put axon models at the specified distances and compute their minimum required stimulation amplitude to elicit an action potential. An example for determining the required stimulation amplitude range between a distance of 2 mm and 4 mm to the active electrode contact is given by
\newpage
\begin{lstlisting}[language=python]
def create_single_threshold_point(distance, axon, dbs_signal, fieldmodel):
	number_of_points = len(axon_min_max.get_points())
	point = Threshold_Point(Threshold_Location(0, 0, number_of_points)
	transform = point.get_transform()
	transform.set_position_vector(np.array([0, distance, 0])
	nodes = transform.coord_transform(axon.get_nodes())
	point.set_coordinate(nodes[0]/2+nodes[-1]/2)
	point.set_axon(axon)
	point.set_dbs_signal(dbs_signal)
	potential = fieldmodel.evaluate_quantity(fieldmodel.get_potential(), nodes)
	point.set_potential(potential)
	return point

ncpus = 1
axon_min_max = Threshold_Volume()
axon_model = Axon({'diameter': 5.7})
point_min = create_single_threshold_point(2e-3, axon_model, dbs_signal, fieldmodel)
axon_min_max.add_point(point_min)
point_max = create_single_threshold_point(4e-3, axon_model, dbs_signal, fieldmodel)
axon_min_max.add_point(point_max)
axon_min_max.compute_optimal_activation(ncpus, 1.0)
stim_min = axon_min_max.get_point(point_min).get_threshold()
stim_max = axon_min_max.get_point(point_max).get_threshold()
\end{lstlisting}

\paragraph{Field Threshold Approximation of the Neural Activation Extent}
Following an approach presented by \href{https://www.ncbi.nlm.nih.gov/pubmed/25350910}{\r{A}str\"om, et al, 2015}, field threshold values of field quantities can provide a good approximation the neural activation extent. Nevertheless, using pre-determined field threshold values for a certain axon diameter and field model can result in deviations if the model parameters of the field model are varies. To overcome these deviations, the field treshold values in FanPy are determined by computing the threshold-distance relationship for the corresponding field and axon model and using this relationship to determine the location at which the axon gets activated by a given stimulation amplitude of the applied DBS signal. The following image shows threshold-distance relationship for axons with varying fiber diameter.
\begin{figure}[!th]
    \centering
		\includegraphics[width=0.8\textwidth]{./pics/threshold_distance_relationship.png}
		\caption{Threshold-distance relationship for axon fibers with varying fiber diameters $f_d$.}
		\label{fig:adaptivealgorithm}
\end{figure}
In FanPy an example code to use the field apprximation approach to estimate the neural activation extent by the iso-volume of an electric potential threshold value is given by

\begin{lstlisting}[language=python]
npcus = 1
neuronparams = {
                'neuronparams': {
                    'g_step': g_step,
                    'rot_degree_step': rot_degree_step,
                    'dbs_signal': dbs_signal,
                    'fiberD': fiberD,
                    'approx_by_field': True,
					'approx_by_field_quantity': 'potential',
                    'stim_min': 2.0,
                    'stim_max': 4.0,
                    'symmetric': symmetric,
                    }		
			}
vta = vtautil.compute_vta(ncpus, neuronparams)
\end{lstlisting}

\end{document}