1 Introduction

One-third of individuals with focal epilepsy continue to have seizures despite optimal medical management. These patients are candidates for resection if the epileptogenic zone (EZ) can be identified. Intracranial depth electrodes may be implanted in the brain to record electroencephalographic (EEG) signals indicative of epileptic activity in both deep and superficial regions of interest (ROIs) within the cortex that have been identified as potential EZ. Implanted electrodes are also used for stimulation studies to map eloquent areas (e.g. motor or sensory cortex) and to determine whether a safe resection may be made that removes the EZ without compromising eloquent cortex.

Electrode trajectories must avoid critical structures (blood vessels, sulci) to minimise the risk of haemorrhage. In current clinical practice each electrode trajectory, defined by a brain target and a skull entry, is planned by visual inspection of preoperative imaging data. Trajectories are placed to attain deep and superficial ROIs, avoid critical structures, sample grey matter (GM), and avoid conflict between electrodes. Maximal sampling of GM is preferred, as epileptic seizures arise in GM rather than white matter. This is a complex, time-consuming task typically requiring 2–3 h per case. Trajectory planning algorithms may reduce planning time, improve safety, and increase GM sampling by optimizing trajectories according to quantitative suitability measures.

Single trajectory planning algorithms require the user to specify a target as a point [2, 6, 8] or manual ROI [5]. Trajectories are discarded if surgically infeasible or an unsafe distance from critical structures. Trajectories are then scored with a risk metric and the most suitable trajectory is selected. These methods are limited when applied to epilepsy as (1) the most suitable target point for deep ROIs may not be obvious, (2) multiple electrodes are implanted and must be placed to avoid conflicts between electrodes and to maximize GM sampling.

Trajectory planning algorithms designed for intracranial depth electrodes allow users to specify target and entry ROIs [4, 7]. de Momi et al. [4] determined trajectories by randomly sampling from user selected target and entry ROIs. Zelmann et al. [7] determined candidate target points within the hippocampus and amygdala by sampling a Gaussian distribution defined on a ROI distance map where targets far from the surface were preferentially sampled. For both methods, the best trajectory in terms of entry angle with the skull, avoidance of critical structures, and no conflicts between electrodes was computed. Zelmann et al. [7] also had a constraint to maximize GM sampling.

We present anatomically-driven multiple trajectory planning (ADMTP) to compute electrode trajectories. ADMTP improves the state-of-the-art by (1) requiring the user to only input anatomical names for deep and superficial cerebral ROIs, (2) selecting target points within deep ROIs according to safety, measured as the distance from critical structures, and (3) computing the best trajectories to attain superficial ROIs, avoid critical structures, improve GM sampling, and avoid conflicts between electrodes.

2 Methodology

ADMTP requires brain parcellation and segmentation of critical structure (Sect. 2.1). From a list of user defined ROIs ADMTP finds an implantation plan \(V( N ) = [ v_{1}, \ldots , v_{N} ]\), where \(v_{n}\) is the trajectory for the n th electrode. For \(v_{n}\), M candidate target points \(T_{n,i} : i \in \{1, \ldots , M\}\) are calculated as described in Sect. 2.2. Next ADMTP determines P entry points \(E_{n,j} : j \in \{1, \ldots , P\}\) for each \(T_{n,i}\) and then a combination of \(T_{n,i}\) and \(E_{n,j}\) for all N electrodes are computed to avoid electrode interference as detailed in Sect. 2.3.

2.1 Regions of Interest and Critical Structure Extraction

Brain parcellation and segmentation of GM, arteries, veins, sulci, and the skull is performed as follows. Veins and arteries are segmented from CT angiography or T1 weighted MRI with gadolinium enhancement using multi-scale, multi-modal tensor voting [9]. The skull is segmented from CT using thresholding. Geodesic Information Flows (GIF) [3] is used to parcelate the brain on T1 weighted MRI, giving 208 ROIs. GM and sulci are obtained from the brain parcellation. Figure 1 shows an implantation plan with deep and superficial ROIs.

Fig. 1.
figure 1

Implantation plan with 10 electrodes. (a) Deep ROIs: amygdala (cyan), hippocampus (yellow), anterior insula (brown), transverse temporal gyrus (blue), posterior (orange) and middle cingulate (purple), posterior (green) and anterior medial orbital gyrus (mauve). (b) Superficial ROIs: middle frontal (blue), superior frontal (purple), middle temporal (light pink), and superior temporal (dark pink) lobes and supramarginal (light orange), angular (light green), and precentral (dark green) gyri.

2.2 Candidate Target Point Selection

For a ROI M candidate target points \(T_{n,i} : i \in \{1, \ldots , M\}\) are computed for the n th electrode. First, a target risk image is defined as \(\mathcal {C}_{t} = [C, f_{t}(c): c \in C]\) where C is the grid of image voxel locations and \(f_{t}(c)\) is the target risk score for the given voxel c. \(f_{t}(c)\) is computed as,

$$\begin{aligned} f_{t}(c) = {\left\{ \begin{array}{ll} 1 ,&{} \text {if } c \notin \varOmega _{roi} \\ 1, &{} \text {if } c \in \varOmega _{cri} \\ w_{roi} * \bigg ( 1 - \dfrac{f_{roi}(c)}{\max (f_{roi})} \bigg ) + w_{cri} * \bigg ( 1 - \dfrac{f_{cri}(c)}{\max (f_{cri})}\bigg ), &{} \text {else} \end{array}\right. } \end{aligned}$$
(1)

\(\varOmega _{roi}\) and \(\varOmega _{cri}\) are the set of voxels in the deep ROI and critical structures, respectively. \(f_{roi}(c)\) is the distance between c and the closest surface point on the deep ROI, calculated using a bounding volume hierarchy (BVH) as in [8]. Similarly, \(f_{cri}(c)\) is the distance between c and the closest surface point on all critical structures. \(w_{roi}\) and \(w_{cri}\) control the relative importance of placing the target within the deep ROI and avoiding critical structures, respectively. Figure 2 displays \(\mathcal {C}_{t}\) for a hippocampus, red corresponds to high values of \(f_{t}(c)\) and green to low values. Target points are selected from \(c \in \mathcal {C}_{t}\) by calculating local minima using the watershed algorithm [1] then sampling the M lowest values of \(f_{t}(c)\) with a distance of at least \(d_{tar}\) between every target point \(T_{n,i}\).

Fig. 2.
figure 2

(a)Axial, (b) sagittal, and (c) coronal views of \(\mathcal {C}_{t}\) (high values of \(f_{t}(c)\) are red, low values are green) for the hippocampus (blue) with blood vessels (cyan) and trajectories determined by ADMTP (purple, pink).

2.3 Automated Trajectory Planning

Entry points, defined as \(E_{n,j} : j \in \{1, \ldots , P\}\), are computed by using all vertices on the skull mesh, roughly 10, 000 points sampled every 0.2 mm3. Potential trajectories, \(\overline{T_{n,i}E_{n,j}} : i \in \{1, \ldots , M\}, j \in \{1, \ldots , P\}\), are then removed from consideration using a modified approach of Zombori et al. [8] as follows:

  1. 1.

    \(\overline{T_{n,i}E_{n,j}}\) is longer than \(d_{len}\).

  2. 2.

    The angle between \(\overline{T_{n,i}E_{n,j}}\) and the skull normal is greater than \(d_{ang}\).

  3. 3.

    \(\overline{T_{n,i}E_{n,j}}\) does not traverse the superficial ROI.

  4. 4.

    \(\overline{T_{n,i}E_{n,j}}\) intersects a critical structure (arteries, veins, or sulci).

After excluding trajectories based on these hard criteria there are typically 1, 000–5, 000 potential trajectories per electrode. The remaining trajectories have a risk score \(R_{n,i,j}\) and GM ratio \(G_{n,i,j}\) calculated. \(R_{n,i,j}\), a measure of the distance to critical structures, is computed as,

$$\begin{aligned} R_{n,i,j} = \dfrac{\displaystyle \int _{E_{n,j}}^{T_{n,i}} d_{risk} - (f_{cri}(x)-d_{safe} ) dx}{(d_{risk}-d_{safe}) * length}, \end{aligned}$$
(2)

where trajectories with \(f_{cri}(x)\) closer than \(d_{safe}\) have the highest risk (\(R_{n,i,j}=1\)) while \(f_{cri}(x)\) farther than \(d_{risk}\) have no risk (\(R_{n,i,j}=0\)).

\(G_{n,i,j}\) measures the proportion of electrode contacts in GM. For each electrode, Q contacts with a recording radius of \(p_{r}\) are spaced at even intervals \(p_{q}\) along the trajectory. \(G_{n,i,j}\) is calculated as,

$$\begin{aligned} G_{n,i,j} = \dfrac{\displaystyle \sum _{q = 1}^{Q} (H[f_{gm}(p_q - p_r) ] + H[f_{gm}(p_q) ] +H[ f_{gm}(p_q + p_r) ] }{3 * Q}, \end{aligned}$$
(3)

where \(f_{gm}(\cdot )\) is the signed distance from the GM surface and \(H[\cdot ]\) is the Heaviside function, with values of 1 inside GM and 0 outside.

Table 1. The following values were set by a consensus of 3 neurosurgeons: the most oblique angle drillable, \(d_{ang}\), the minimum safe distance from blood vessels, \(d_{safe}\), the distance at which there is no risk, \(d_{risk}\), the minimum distance between electrodes, \(d_{traj}\). A commonly used electrode configuration determined: electrode length, \(d_{len}\), the number of contacts, Q, the interval between contacts, \(p_q\), the contact sample radius, \(p_r\). The following values were set empirically: the number of candidate targets, M, the minimum distance between candidate targets, \(d_{tar}\), and the relative importance of sampling the ROI, \(w_{roi}\), and avoiding critical structures, \(w_{cri}\).

Each trajectory is assigned a weighted score \(S_{n,i,j}\) computed as \(S_{n,i,j} = 10 * R_{n,i,j} + G_{n,i,j}\), where 10 was determined empirically so low risk is prioritized over a high GM ratio.

The final implantation plan V(N) is found by optimizing,

$$\begin{aligned} S_{total}&= \mathop {\hbox {argmin}}\limits _{V(N)} {\bigg ( \dfrac{1}{N} \sum _{n=1}^{N} S_{n,i,j} \bigg )} \nonumber \\&\text { s.t. } D(\overline{T_{n,i} E_{n,j}}, \overline{T_{k,i} E_{k,j}} ) > d_{traj} : \forall n, \forall k \in \{ 1, \ldots , N\}, n \ne k. \end{aligned}$$
(4)

where \(d_{traj}\) specifies the minimum distance between trajectories that do not conflict. Due to the constraint \(d_{traj}\) if the user selects multiple electrodes with the same ROI ADMTP will find unique targets. For an implantation plan there are typically 7–12 electrodes, each with 1, 000–5, 000 potential trajectories, representing approximately \(1 \times 10^{21}\) possible combinations, hence, a depth-first graph search strategy is used to calculate a feasible implantation plan. If no combination of trajectories exists which satisfies \(d_{traj}\) ADMTP returns the plan with the largest distance between trajectories.

Fig. 3.
figure 3

Measures of suitability for trajectories determined by manual planning (plotted on the X axis) versus ADMTP (plotted on the Y axis) are shown for (a) angle with re-spect to the skull surface normal, (b) risk score, (c) distance to the nearest critical structure, and (d) GM ratio. Red points are the center of mass for each measure. For (a) angle and (b) risk score points below the diagonal correspond to ADMTP giving the preferred result. For (c) distance to critical structures and (d) GM ratio points above the di-agonal correspond to ADMTPgiving the preferred result.

3 Experimental Design and Results

3.1 Experimental Design

ADMTP was evaluated on retrospective data from 20 patients with refractory focal epilepsy who underwent intracerebral electrode implantation. Each patient had 7–12 electrodes, 186 in total. Manual plans were determined by the consensus of two neurosurgeons. Deep and superficial ROIs were identified by the same neurosurgeons and ADMTP was evaluated using the parameters in Table 1.

3.2 Trajectory Suitability

Trajectories were assessed by angle with respect to the skull surface normal, risk score, distance to nearest critical structure, and GM ratio. In Fig. 3 each point corresponds to one trajectory with the manual plan value plotted on the X axis and the ADMTP value plotted on the Y axis. The red point represents the center of mass for each measure. Points below the diagonal have a lower value for ADMTP compared to manual plans which represent ADMTP giving the preferred value for angle and risk score; points above the diagonal represent ADMTP giving the preferred value for critical structures distance and GM ratio. A two-tailed Student’s t-test evaluated the statistical significance between values determined by ADMTP and manual plans where the null hypothesis was that the methods return similar values.

ADMTP found a more feasible entry angle in 96 / 186 trajectories (\(p < 0.01\)) and increased GM sampling in 104 / 186 trajectories (\(p > 0.01\)). ADMTP found trajectories that were safer, in terms of reduced risk score and increased distance to the closest critical structure in 145 / 186 trajectories (\(p < 0.01\)).

Fig. 4.
figure 4

Manual (pink) and ADMTP (blue) trajectories are shown with veins (cyan), skull (opaque white), and with (a) the cortex (peach) and (b) no cortex.

Table 2. Computation time for each step in ADMTP reported in seconds.

3.3 Implantation Plan Suitability

Suitability of implantation plans were assessed by (1) distance between trajectories and (2) the ratio of unique gyri sampled to total number of electrodes. Ideally each electrode samples a unique gyrus, corresponding to a ratio of 1. ADMTP has an median distance between trajectories of 35.5 mm (5.2–124.2 mm) compared to manual plans with an median of 34.2 mm (1.3–117.5 mm). Manual plans and ADMTP have 12 trajectory pairs separated by less than \(d_{traj}\) (<10 mm). Manual plans had an average gyri to electrode ratio of 0.96 (0.86-1) while ADMTP has a ratio of 0.92 (0.75–1). ADMTP has larger distances between trajectories but manual planning still provides preferred coverage in terms of more unique gyri sampled and consistent coverage of the cortex. Figure 4 displays a manual (pink) and ADMTP (blue) plan.Figure 4(b) shows target point differences within the cortex, and Fig. 4(a) shows entry point differences. ADMTP often finds trajectories within the same gyrus as manual planning.

3.4 Computational Efficiency

The computational efficiency of ADMTP was assessed for computing candidate target points, trajectories, and ADTMP (combination of both steps). Computations were performed on a computer with a Intel(R) Xeon(R) 12 core CPU 2.10 GHz with 64.0 GB RAM and a single NVIDIA Quadro K4000 4 GB GPU. Table 2 reports computation time. ADMTP took 61.14 s (15.43–279.20 s) with most of the computation time spent on trajectory planning.

4 Concluding Remarks

We presented an anatomically driven multiple trajectory planning (ADMTP) algorithm for calculating intracerebral electrode trajectories from anatomical regions of interest (ROIs). Compared to manual planning, ADMTP lowered risk in \(78\,\%\) of trajectories and increased GM sampling in \(56\,\%\) of trajectories. ADMTP was evaluated on quantitative measures of suitability, however, a qualitative analysis is necessary to assess the clinical suitability of ADMTP. Future work is required to ensure ADMTP provides trajectories that sample unique gyri. ADMTP efficiently calculates (>5 min) safe trajectories.