Automated multiple trajectory planning algorithm for the placement of stereo-electroencephalography (SEEG) electrodes in epilepsy treatment

Purpose About one-third of individuals with focal epilepsy continue to have seizures despite optimal medical management. These patients are potentially curable with neurosurgery if the epileptogenic zone (EZ) can be identified and resected. Stereo-electroencephalography (SEEG) to record epileptic activity with intracranial depth electrodes may be required to identify the EZ. Each SEEG electrode trajectory, the path between the entry on the skull and the cerebral target, must be planned carefully to avoid trauma to blood vessels and conflicts between electrodes. In current clinical practice trajectories are determined manually, typically taking 2–3 h per patient (15 min per electrode). Manual planning (MP) aims to achieve an implantation plan with good coverage of the putative EZ, an optimal spatial resolution, and 3D distribution of electrodes. Computer-assisted planning tools can reduce planning time by quantifying trajectory suitability. Methods We present an automated multiple trajectory planning (MTP) algorithm to compute implantation plans. MTP uses dynamic programming to determine a set of plans. From this set a depth-first search algorithm finds a suitable plan. We compared our MTP algorithm to (a) MP and (b) an automated single trajectory planning (STP) algorithm on 18 patient plans containing 165 electrodes. Results MTP changed all 165 trajectories compared to MP. Changes resulted in lower risk (122), increased grey matter sampling (99), shorter length (92), and surgically preferred entry angles (113). MTP changed 42 % (69/165) trajectories compared to STP. Every plan had between 1 to 8 (median 3.5) trajectories changed to resolve electrode conflicts, resulting in surgically preferred plans. Conclusion MTP is computationally efficient, determining implantation plans containing 7–12 electrodes within 1 min, compared to 2–3 h for MP.


Introduction
Between 20 and 40 % of focal epilepsy patients are refractory to antiepileptic medications [13]. Such patients are candidates for curative surgery, which aims to resect the epileptogenic zone (EZ) that generates seizures [3]. In about 25 % of surgical candidates, the EZ cannot be inferred from noninvasive imaging data, and intracranial electroencephalography (EEG) is needed to identify the EZ [7].
Stereo-EEG (SEEG) records EEG signals via depth electrodes surgically implanted in the brain. SEEG electrodes record from a 1-cm core around the cerebral entry to the distal end (target) that may be placed in hippocampus, amygdala, or midline or neo-cortex in temporal, frontal, parietal, or occipital lobes. Electrode implantation carries a risk of haemorrhage, neurologic deficit, and infection [4].
Preoperative planning of electrode trajectories, defined by the target and the skull entry point, can minimise implantation risk by ensuring electrodes avoid critical structures (e.g. arteries, veins, sulci) and conflicts between electrodes. Plan-ning may also improve the efficiency of the SEEG recording by ensuring electrodes pass through the maximal amount of grey matter (GM), GM being the component of brain tissue that generates seizures. Current clinical practice for planning electrode trajectories involves manual evaluation of trajectories in series. This is a complex, time-consuming task requiring: 1. Integrating information across imaging modalities to locate critical structures, GM, and targets. 2. Optimising several criteria for each trajectory to sample the target, avoid critical structures, and obtain a suitable angle to traverse the skull. 3. Adjusting trajectories to maximise GM capture and avoid electrode conflicts, whereby two electrodes may contact each other. Placing a new electrode may require adjusting previously planned trajectories.
Computer-assisted planning algorithms can reduce planning time by calculating quantitative measures of trajectory suitability. These measures can be used to select the best trajectory (automated planning) or inform manual trajectory selection (assisted planning).
We present an automated multiple trajectory planning (MTP) algorithm that calculates a combination of trajectories, or plan, for a set of targets. Trajectories are assessed on proximity to critical structures (risk score) and sampling of GM [GM-white matter (GM-WM) ratio]. Our MTP algorithm uses dynamic programming to reduce the search space and a depth-first search to find a plan whereby each electrode trajectory is surgically feasible, does not interfere with other trajectories, avoids critical structures, and maximises GM sampling. MTP is integrated into the EpiNav TM software platform [29] to enable manual trajectory assessment.
The remainder of the manuscript is organised as follows. The third section describes the previous work in computerassisted trajectory planning. The fourth section describes our MTP algorithm. The fifth section describes the evaluation of our MTP algorithm. The sixth section discusses MTP, and the seventh section provides concluding remarks.
Assisted planning methods aim to reduce the time and complexity of manual trajectory selection by displaying mea-sures of risk for potential trajectories [11,16,20,22]. [16] displayed a heat map corresponding to the minimum distance to critical structures for potential entry points. Similarly, [11] displayed an entry point safety map, safety being related to distance from critical structures. [20] reduced computation time for safety maps using graphical processing units (GPUs) to enable real-time user interaction. [22] displayed a cumulative risk, the summation of distance from critical structures along the trajectory.
Single trajectory planning algorithms automatically determine the best trajectory for one electrode given a specific quantitative measure. [25] assessed trajectories using a risk score calculated by summing traversal costs, where regions to be avoided had a high traversal cost, along the trajectory. Similarly, [10] summed traversal costs along the trajectory but added a penalty for trajectories near blood vessels to reduce the risk of haemorrhage. [1] assessed trajectories by first removing potential trajectories that were an unsafe distance from critical structures. The remaining trajectories were assessed by a weighted sum of (1) the minimum distance to critical structures and (2) the cumulative distance from all critical structures. [23] calculated a traversal cost by first computing a per pixel risk score, based on distance to critical structures. They then determined two risk scores: (1) maximum risk along the trajectory and (2) a summation of the risk along the trajectory. The user could select which of these two risk scores to use. [8] developed a generic optimisation algorithm for trajectory planning, allowing a user to define a set of hard constraints, rules that must not be violated, and soft constraints, rules that could be minimised. The generic optimisation eliminates trajectories that violate the hard constraints and then finds the trajectory which minimises the summation of soft constraints. [2] used a similar approach, defining hard constraints, specific entry points and avoiding critical structures (ventricles, blood vessels, sulci), and soft constraints, minimising overlap with the caudate and GM. [14] combined 6 soft constraints and 2 hard constraints to define a weighted cost function that determined the best trajectory for targeting the subthalamic nucleus. [26] developed an algorithm that maximises distances from critical structures and GM sampling. [24] developed an algorithm that optimises a weighted sum of the trajectory distance to critical structures so that blood vessels, with a high weight, are always avoided, while WM tracts and regions of cortical function, with a low weight, may be traversed if no other path exists. Constraints in [24] are specific to placing electrode for DBS and may not be generalisable to SEEG electrode implantation.
Multiple trajectory planning algorithms determine the best combination of trajectories, or plan, for multiple electrodes. Multiple trajectory planning not only takes into account the quantitative measures for individual electrodes but also that electrodes must not contact each other. [26] optimised three electrodes with targets in the amygdala, anterior and posterior hippocampus so that risk was minimised and there were no conflicts between electrodes. In [26] potential entry points were constrained by the use of target specific entry map priors. A more extensive evaluation of this method on 37 patients was presented in [27] in which increased GM sampling by optimising the number of electrode contacts in GM was demonstrated.
[5] presented a multiple trajectory planning algorithm that operates serially. The algorithm sets the first electrode to the best trajectory, and additional electrodes are set by removing trajectories that interfere with existing electrodes, and then selecting the best trajectory. This algorithm is dependent on trajectory order and may not find an optimal plan. [6] overcame these limitations by evaluating all potential plans and returning the best plan. However, enumerating all potential plans is computationally expensive when evaluating many targets or many potential trajectories per target. To reduce the number of potential plans, [6] reduced potential entry points by randomly sampling within a user-specified region.
Previous work from our group presented an automated entry point search algorithm [29] in which all points on the skull were potential entry points, removing the need to manually specify entry regions or define entry map priors for each target. All points on the skull result in 2000-10,000 potential entry points, and for eight targets this corresponds to a minimum of 1E26 potential plans.
We present a MTP algorithm that combines dynamic programming to reduce the number of potential plans, and depth-first searching, to find a suitable plan. Trajectories are assessed by risk score, measured as the cumulative distance to blood vessels from the trajectory [29], and GM sampling, measured as the proportion of electrode contacts that are in GM. Our algorithm is an improvement over the current state of the art in that it: (1) finds a combination of electrode trajectories with no limitations on the number or order of electrodes and (2) is computationally efficient, finding a plan with 7 to 12 electrode in under one minute, and (3) is integrated in the EpiNav TM software platform to enable manual assessment.

Multiple trajectory planning algorithm
A trajectory is defined as v = {T, E, R, G} where T is the target in the brain, E is the entry on the skull, R is the risk score, and G is the GM-WM ratio. For a set of N targets a plan is defined as V (N ) = {v 1,a 1 , . . . , v N ,a N } : a i ∈ {1, . . . , M i }, i ∈ {1, . . . , N } where M i is the number of potential trajectories for the i th target. The plan V (N ) is defined such that each trajectory attains one of N targets. MTP finds a plan V min (N ) that attains all targets, minimises R, maximises G, and avoids conflicts between electrodes.
Prior to trajectory planning, segmentation of the skull and critical structures is performed as described in the section "Critical structure extraction". For each target T i , a risk score R i,a i , that quantifies proximity to blood vessels, and a GM-WM ratio G i,a i , that quantifies GM sampling, are calculated as described in the section "Single trajectory planning algorithm". The MTP algorithm described in the section "Multiple trajectory planning algorithm" calculates V min (N ). The EpiNav TM software platform enables manual assessment of V min (N ) as described in the section "Plan visualisation and assessment".

Critical structure extraction
Automated trajectory planning is dependent on accurately segmenting critical structures (arteries, veins, and sulci), GM, and the skull surface. Algorithms chosen for these tasks are currently being used in the clinic to generate 3D models for manual trajectory planning and have been used in presurgical patient evaluation for over 4 years [19,21]. Blood vessels are segmented with a customised vessel extraction tool [30] from CT angiography, 3D phase contrast MRI, or T1-weighted MRI with gadolinium enhancement. GM and the cortex were segmented using FreeSurfer [9]. Sulci were extracted from the cortex surface. Figure 1c illustrates an example segmentation for veins (cyan), arteries (red), and sulci (peach).
Skull segmentation with template registration constrains entry points to regions suitable for implantation. A patientspecific skull is segmented from a CT scan using thresholding Fig. 1 Example study displaying a skull segmentation (white), b skull template (yellow) that excludes surgically infeasible regions, and c veins (cyan), arteries (red), and sulci (peach) with skull (semi-transparent white) and skull template (semi-transparent yellow) and morphologic dilation to ensure a fully connected surface. The template skull is aligned to the patient skull using Iterative Closest Points (ICP) [28] to minimise the distance between the two surfaces. The template skull excludes regions inappropriate for implantation such as the face, ears, and regions inferior to the transverse sinus. Figure 1 shows an example patient (white) and template (yellow) skull.

Single trajectory planning algorithm
Previous work from our group [29] presented real-time automated single trajectory planning (STP) for a target T i . STP calculates potential entry pointsÊ i,a i : a i ∈ {1, . . . , M i } by considering trajectory length and entry angle (described in the section "Entry point search"). Next trajectories, defined asÊ i,a i T i , that intersect critical structures (arteries, veins, or sulci) are removed from consideration. Then at evenly spaced points along each trajectory x ∈ E i,a i T i the distance to the nearest blood vessel f crit (x) is found (described the section "Bounding volume hierarchy (BVH) for trajectory evaluation"). Finally for E i,a i T i a risk score R i,a i , computed from f crit (x), and a GM-WM ratio G i,a i are calculated (described in the section "Trajectory ranking"). A stratified ranking algorithm sorts trajectories by first minimising R i,a i and then maximising G i,a i .

Entry point search
Potential entry pointsÊ i,a i : a i ∈ {1, . . . , M i } are identified by considering all vertices in the template skull mesh. E i,a i are removed from consideration based on the following criteria: 1. Trajectory length The length ofÊ i,a i T i must be shorter than d length , the maximum electrode length.
2. Entry angle The angle betweenÊ i,a i T i and the skull normal must be less than d angle , the angle that can be accurately drilled.
Calculating these exclusion criteria forÊ i,a i T i is computationally inexpensive; hence, it is practical to removeÊ i,a i that do not meet these criteria first.

Bounding volume hierarchy (BVH) for trajectory evaluation
Each trajectoryÊ i,a i T i is tested for intersection with critical structures (arteries, veins, sulci) using a bounding volume hierarchy (BVH) to enable real-time calculation. Trajectories that intersect these structures are removed from consideration. All remaining trajectories are sampled at 128 evenly spaced points x such that x ∈ E i,a i T i . For every x, the distance to the nearest blood vessel (arteries, veins) f crit (x) is calculated. BVH construction and traversal are described below.
Bounding volume hierarchy construction For each critical structure (arteries, veins, and sulci) a BVH is constructed as in [12]. Each triangle in the surface is assigned a 30-bit Morton code [15], calculated by combining the 10-bit Morton code of each triangle vertex coordinate. An efficient bit-wise sorting of the triangles is performed using the Morton codes. The BVH is created by iteratively splitting triangles according to the highest different bit between Morton codes. This is repeated until each leaf node contains one triangle. Finally, for every node a bounding box is calculated. For each leaf node the bounding box is calculated as the smallest rectangle that contains the triangle. The bounding box for all other nodes is the union of the bounding boxes of their children nodes.
Bounding volume hierarchy traversal BVH traversal detects collision ofÊ i,a i T i with each critical structure. Initially the top BVH node is added to the queue. IfÊ i,a i T i intersects the bounding box of the first node in the queue, its children nodes are added to the queue. Once a leaf node is reachedÊ i,a i T i is removed from consideration if it intersects the triangle of the leaf node. For the remaining trajectories, the closest distance between each point x ∈ E i,a i T i and the jth critical structure f j (x) is calculated. For this computation sulci are not included. Initially the top BVH node is added to the queue and f j (x) = ∞. The first node in the queue is removed, and the distance between each child node and its bounding box f bb (x) is calculated. For a point inside the bounding box f bb (x) = 0. If f bb (x) < f j (x) the node is added to the queue so that the first node corresponds to the smallest value of f bb (x). For a leaf node the distance between x and the tri- This is repeated until no nodes are in the queue. After all critical structures have been evaluated, the closest distance is calculated as f crit (x) = arg min j ( f j (x)).

Trajectory ranking
Entry points that meet all hard constraints, E i,a i : a i ∈ {1, . . . , M i }, are ranked by risk score R i,a i , a measure of cumulative distance from blood vessels, and GM-WM ratio G i,a i , a measure of GM capture.
Risk score The risk score R i,a i measures cumulative distance to blood vessels. The trajectory E i,a i T i has a high risk if the nearest critical structure is less than a "Safety Margin", determined by the user-defined value d safety . If E i,a i T i has a . a Initial set of potential entry pointsÊ i,a i obtained from the skull template (semi-transparent white). b Opaque white patch corresponds toÊ i,a i after excluding trajectories due to length or angle. c Opaque white patch corresponds toÊ i,a i after excluding trajectories due to critical structure intersection. d Coloured patch corresponds to risk scores R i,a i with low (0-green) to high (1-red) risk distance to the nearest critical structure greater than a "Risk Zone", determined by the user-defined value d risk , it has no potential risk.
The cumulative distance of risk along E i,a i T i is calculated as, where f crit (x) is the distance between x and the nearest critical structure. For normalisation purposes if f crit (x) > d risk , then f crit (x) = d risk so that the final value the x contributes to the risk score is zero. If f crit (x) < d safety , then automatically R i,a i = 1, representing the highest risk. The final risk R i,a i is normalised to the range [0, 1], where 0 corresponds to no risk and 1 corresponds to the highest risk. R i,a i is calculated as, where length is the length of E i,a i T i . Figure 2d displays R i,a i as a heat map from low (0-green) to high (1-red) risk.
Grey matter-white matter ratio GM-WM ratio measures the proportion of electrode contacts in GM. GM-WM ratio corresponds to the SEEG efficiency for each trajectory as GM generates seizures. For each trajectory E i,a i T i a set of J contact points, c j : j ∈ {1, . . . , J } each with a sampling radius c r are defined. Each contact point is assessed if c j ± c r is located in GM. GM-WM ratio is calculated as, Stratified ranking Trajectories are first ranked by R i,a i so that v i,1 has the lowest risk. Next trajectories are placed into K histogram bins so the k th bin contains v i,a i : where D(·, ·) is the minimum Euclidean distance between two line segments. A depth-first search with dynamic programming to limit potential plans is used to calculate V min (N ) as described in the sections "Depth-first search algorithm" and "Dynamic programming for determining potential combinations".

Depth-first search algorithm
Algorithm 1 iteratively (1) calculatesV min (n), a suitable plan for n trajectories and (2) rejectsV min (n) if electrodes Algorithm 1: Depth-first search to find the optimal plan for N targets.
conflict. At each iteration of Algorithm 1 a set of low risk plans, defined as V p (n) = [V 1 (n), . . . , V q (n)], is calculated. For the lowest risk planV min (n) ∈ V p (n) the function D all (·) detects conflict between electrodes by finding the minimum Euclidean distance between all pairs of trajectories inV min (n) (i.e. min(D(E i,a i T i , E j,a j T j ) : ∀i, ∀ j ∈ {1, . . . , n}, i = j)). When an electrode conflict is detected V p (n) is updated as described in the section "Dynamic programming for determining potential combinations". Oncê V min (n) has no conflicts, n ←− n + 1 and the algorithm continues until V min (N ) is found.

Dynamic programming for determining potential combinations
For each target T i , trajectories are ranked so v i,1 is the best trajectory (in the section "Trajectory ranking"). Initially, V p (n) ←− [V 1 ] where V 1 = {v 1,1 } which corresponds to adding a potential plan containing the best trajectory for the 1st electrode. At the next step, n = 2, V p (n) ←− [{V min (n − 1), v n,1 }], which corresponds to adding the best trajectory for the 2nd electrode. Algorithm 1 proceeds in this manner until electrode conflict is detected.
Algorithm 1 continues in this manner until V min (N ) is found.

Plan visualisation and assessment
Clinicians must be able to visualise and assess trajectory feasibility. We have developed the EpiNav TM software platform to aid manual assessment with: (1) quantitative measures of trajectory suitability, (2) a trajectory profile, and (3) a probe eye view. Figure 3 displays an example layout of the EpiNav TM software platform.

Quantitative measures
Four measures of trajectory suitability, length, angle, risk score, and GM-WM ratio, are displayed in the EpiNav TM platform. Trajectory length is calculated as the length of E T , where E is the entry point and T is the target. Trajectory angle is calculated as the angle of E T with respect to the skull normal. Trajectory risk score is calculated as in Equation 2, and GM-WM ratio is calculated as in Equation 3.

Trajectory profile
The trajectory profile provides a (a) risk profile and (b) GM profile. The risk profile displays f crit (x) (described in the section "Trajectory ranking") for x ∈ E T . The colour of f crit (x) corresponds to the closest blood vessel at x. The GM profile displays the colour corresponding to the tissue type (GM or WM) each electrode contact is located in. Regions of the trajectory outside the cortex are shown in the electrode colour. Figure 3 displays the risk profile in the upper right panel. Red corresponds to trajectory regions closest to an artery, cyan to trajectory regions closest to a vein. The red line indicates the "Safety Margin" (d safety ). The black line indicates the position of the probe eye view (Fig. 3, bottom right panel), black text shows the distance between the probe eye view and T , and red text shows the distance to the nearest critical structure ( f crit (x)).

Probe eye view
The probe eye view displays an image perpendicular to the trajectory. This allows the user to navigate along the trajectory and assess proximity to critical structures. The probe eye view is generated by finding the geometric plane perpendicular to x ∈ E T . Nearest neighbour interpolation calculates the intensity value for each pixel in this plane. Figure 3 displays a probe eye view in the bottom right panel.

Dataset description
Evaluations were performed on retrospective data from 18 patients with medically refractory epilepsy who underwent SEEG implantation. All patients had unilateral implantations with between 7 and 12 electrodes for a total of 165 electrodes. All studies involving human participants were in accordance with the ethical standards of the institutional and/or national research committee and with the 1964 Helsinki Declaration and its later amendments or comparable ethical standards. For this type of study formal consent is not required.

Experimental design
Each electrode target was manually determined by an expert neurosurgeon relying on conventional 2D visualisation. Targets remained fixed across planning methods.
Parameters for STP and MTP were set as described in Table 1. These parameters were determined according to a panel of three expert neurosurgeons. d length , d angle , and d traj are set according to surgically feasible approaches and available electrodes. d safety is set to 3.0 mm corresponding to the minimum accuracy achievable with the neuronavigation system [18] as critical structures closer than 3.0 mm may be compromised. d risk is set to 10 mm, a distance at which there is no potential to compromise critical structures. J , c j , and c r were determined by choosing a common SEEG electrode configuration.
Trajectories were assessed by length, angle, GM-WM ratio, and risk score as described in the section "Quantitative measures". Distance to the closest sulci was also calculated in a similar fashion as distance to blood vessels described in the section "Bounding volume hierarchy (BVH) for trajectory evaluation". These measures of trajectory suitability are biased in that STP and MTP use these measures to calculate electrode trajectories. Although the quantitative measures are biased, they assess whether a method is able to optimise multiple surgical constraints simultaneously. A qualitative analysis by an expert neurosurgeon who was blinded to the plan origin was used to complement the quantitative analy-

Experiment 1: trajectory ranking strategies
We compared trajectories determined by risk score [29] and stratified ranking described in the section "Trajectory ranking". Risk score ranking finds the trajectory with the lowest risk score, while stratified ranking finds a trajectory that minimises risk score while maximising GM-WM ratio. Stratified ranking is expected to increase GM-WM ratio and risk score compared to risk score ranking. A two-tailed Student's t test was used to compare trajectories between the two ranking Fig. 4 a Risk score and b GM-WM ratio for trajectories calculate by risk score ranking (plotted on the Xaxis) versus stratified ranking (plotted on the Y axis). Points above the diagonal represent trajectories where stratified ranking increased the value compared to risk score ranking. For a risk score lower values (points below the diagonal) are preferred. For b GM-WM ratio higher values (points above the diagonal) are preferred methods, with the null hypothesis being that the methods calculate similar trajectories. Figure 4 displays GM-WM ratio and risk score for trajectories calculated with the two ranking methods. Both methods produced trajectories with similar risk scores, and stratified ranking increased the risk score on average 0.02(0.0 − 0.13) and was not statistically significant ( p = 0.497). Stratified ranking increased the GM-WM ratio by a mean of 0.08(0.0 − 0.57) and in 22/165 trajectories by over 0.2. The difference in GM-WM ratio is statistically significant ( p = 5.5×10 −7 ). Stratified ranking improved GM-WM ratio with an insignificant increase in risk score.

Experiment 2: target order independence
We evaluated the effect of target order on MTP. For each plan, target order was randomly selected 5 times and the trajectories for each target were compared. The order targets were considered did not change the final plan, and in all cases the same trajectory was returned.

Experiment 3: planning strategies
We compared our MTP algorithm with (a) manual planning (MP) by an expert neurosurgeon and (b) the STP algorithm described in the section "Single trajectory planning algorithm"

Quantitative assessment
A two-tailed Student's t test was used to compare trajectory measures between MTP and the other methods (MP, STP) with the null hypothesis being the methods return similar trajectories. To account for multiple comparisons (n = 2), a Bonferroni correction is applied; hence, a statistically significant value is α = 0.05/n = 0.025.
All 165 trajectories changed between MP and MTP. , and reduced the distance to the closest sulci in 70/165 trajectories( p = 0.50). For 7 trajectories MTP returned trajectories with a risk score of 1, while MP returned trajectories with a risk score <1; however, these MP trajectories violated the angle constraint (angle > d angle ) (Fig. 6).
Between MTP and STP 69/165 trajectories, these changes were not statistically significant for any measure ( p between 0.03 and 0.60). Although changes in quantitative measure of risk were not statistically significant, there were signifi- cant practical differences in that STP provided implantation plans that could not be surgically implemented due to electrodes being placed too close to each other. In contrast MTP found implantation plans in which each electrode trajectory was placed so that it could be implemented during surgery. MTP calculated a higher risk score than STP in 38/69 trajectories. Additionally, GM-WM ratio decreased in 42/69 trajectories. For one trajectory MTP calculated a much higher trajectory risk score (1 compared to 0.28); however, no low risk trajectories were able to avoid conflicts with other electrodes.

Clinical assessment
A clinical assessment of plan feasibility was performed by a single neurosurgeon, blinded to plan origin. Plans were assessed using EpiNav TM for: 1. Avascularity: each trajectory was assessed with the probe eye view to confirm the absence of nearby blood vessels. Each plan was scored as the ratio of safe, avascular trajectories to all trajectories. 2. Conflicts: each plan was assessed in the volumetric view for conflicts between electrodes, due to either contact or  Table 2 reports the plan feasibility measures. All three methods were effective at finding avascular trajectories for individual trajectories as determined by a neurosurgeon. This Table 2 Measures of plan feasibility for MP, STP, and MTP obtained by one neurosurgeon blinded to plan origin 10 8

Experiment 4: computational time
Computational efficiency of STP and MTP was evaluated. To enable a direct comparison between STP, which calculates the best trajectory for one target, and MTP, which calculates the best trajectories for N targets, the total time to determine all N trajectories was recorded. For STP the computation time is a summation of computation time for each target.
Calculations were performed on a computer with a Intel(R) Xeon(R) 12 core CPU 2.10 with 64.0 GB RAM and a single NVIDIA Quadro K4000 4 GB GPU. Table 3 reports plan computation time. All plans were computed in less than 1 min, and in a clinical setting this will enable the user to make manual adjustments to parameters and trajectories when necessary. Longer computation times were observed for plans with more electrodes (Plans 11 and 13) or electrodes that were placed in close proximity (Plan 7).
Computation time for the preprocessing steps was recorded. GM and cortex segmentation took ≈20 h, surface extraction for the blood vessels and sulci took between 150 and 180 s per structure, and skull segmentation and template registration took between 210 and 260 s. For SEEG electrode implantation patient scans are typically acquired at least 1 week prior to implantation planning; hence, preprocessing steps are not as time sensitive as MTP.

Discussion
Our MTP algorithm is computational efficient, using dynamic programming to consider low risk plans in conjunction with a depth-first search algorithm to find a suitable plan. All plans evaluated containing between 7 and 12 electrodes were calculated in under a minute. Our MTP algorithm resolved electrode conflicts providing more feasible plans compared to STP. Figure 7 displays two plans determined by MP, STP, and MTP. Figure 7c displays a plan in which STP had two electrode conflicts (yellow-blue and pink-purple conflict), and such conflicts prevent the plan being surgically implemented. MTP had no conflicts as shown in Figure 7e. Figure 7d, f illustrates STP and MTP for a plan where one electrode was changed to resolve one electrode conflict (yellow-purple).
Our MTP algorithm has several differences from previously reported multiple trajectory planning algorithms [5,6,26,27]. In terms of target selection, [6] required the user to select a target region from which potential target points were drawn from. [27] constrained target selection to three anatomic regions, amygdala, anterior, and posterior hippocampus. Our MTP algorithm requires the user to specify the target point.
For entry point selection [6] required the user to select an entry region on the skull. [27] defined entry map priors for each anatomic target to constrain potential entry points to surgically feasible regions. Our MTP algorithm uses a generic skull template to constrain entry points making it more flex-ible in the types of electrode trajectories proposed compared to other multiple trajectory planning algorithms. Constraining entry points may be desirable for some electrodes, for example to sample a specific superficial gyrus; however, it may be overly restrictive and result in nonoptimal trajectories for other electrodes.
When calculating trajectories both [6] and [27] sample target and entry regions to obtain a fixed number of trajectory combinations that are then evaluated in terms of risk and electrode conflicts. In contrast, our MTP algorithm considers all possible entry points when determining the trajectories.
Critical structures used to compute the risk score R i,a i (in the section "Critical structure extraction") in this work were blood vessels (arteries and veins). Trajectories that intersected sulci were rejected, but sulci were not included in calculating R i,a i . When analysing MP trajectories, it was found that while blood vessels were avoided by at least 3 mm (median value 5.11 mm), sulci were often much closer (median value 1.57 mm). Based on these results it was determined that maximising distance to sulci was not as important criteria as maximising distance to blood vessels. However, the algorithm presented to compute R i,a i is generalisable to other structures, such as sulci or ventricles, without significantly changing task complexity or expected results. Several and e, f MTP. Each plan shows the skull template (semi-transparent white), critical structures (arteries in red, veins in cyan, and sulci in peach). and trajectories [different coloured entry (arrow) and target (sphere)]. In a, c, and e sulci are not shown, so the electrode configuration can be appreciated other state-of-the-art methods have incorporated sulci avoidance by either taking into account the trajectory angle with respect to the cortex as in [6] or including sulci in the risk metric [27].
The introduction of a skull template (in the section "Critical structure extraction") allows for potential entry points to be limited to surgically feasible regions. The face, ears, and base of the skull are avoided for safety and cosmetic reasons. Figure 8 provides an example where the use of the skull template is necessary to obtain a surgically feasible plan; without the skull template an electrode (orange) would have traversed the posterior fossa inferior to the transverse sinus and penetrated the tentorium cerebelli. However, the skull template as currently implemented is limited. Due to the variability in the position and size of the ears and forehead ICP registration does not always match nonfeasible regions between the patient and template skulls. Individually tailored skull templates would reduce the number of nonfeasible entry points but would increase planning time. Even with individually tailored skull templates our MTP algorithm would still suggest some nonfeasible trajectories. This is due to certain targets having accepted trajectories that neurosurgeons are reluctant to deviate from, even when these trajectories have a higher risk score. In future work, our MTP algorithm will be modified to enable the user to restrict specific electrodes to surgically preferred regions, thereby reducing infeasible trajectories.
In a clinical setting, individual electrode trajectories that are not feasible can be manually adjusted to an appropriate trajectory. However, depending on the location of trajectories, manually adjusting one trajectory may require adjustment of other trajectories to resolve conflicts between electrodes. The EpiNav TM software platform allows the user to fix individual electrode trajectories and rerun MTP on the remaining electrodes to find a suitable implantation plan as described in [17]. Multiple runs of MTP may be required to obtain an implantation plan in which all trajectories are safe and surgically feasible.
A thorough evaluation of our automated multiple trajectory planning algorithm for clinical use has been presented in [17]. In this study, three neurosurgeons compared 18 plans determined by MP and our MTP algorithm. All 18 plans were found to be feasible for clinical implementation. Individual trajectories were found to be safe for clinical implementation in 77.1 %(128/166) of electrodes. 10 %(18/166) of trajectories were found to be unsafe due to incorrect critical structure segmentation. To improve this performance, a more accurate vessel extraction algorithm is necessary to find avascular trajectories. 7 %(12/166) of trajectories were deemed unsafe due to proximity to sulci or the midline, highlighting the need to incorporate sulci for a clinically realistic MTP algorithm.
Stratified ranking allows for a low risk score to be prioritised with GM-WM ratio taken into consideration, provided the risk score does not substantially increase. This additional constraint is important as the goal of electrode implantation is to record EEG signals from GM, which is the site of seizure generation. Currently, the GM-WM ratio is calculated for specific contact points on the electrode for a single-electrode configuration. However, there are over a dozen different configuration of contacts on SEEG electrodes that may be implanted. Future work will include specifying the contact configuration for specific electrodes.

Concluding remarks
We present an automated multiple trajectory planning (MTP) algorithm using depth-first searching with dynamic programming. Our algorithm was evaluated with 18 plans with between 7 and 12 electrodes. Calculation of an implantation plan took on average 19.26(3.86 − 57.38) s. Implantation plans had a lower risk for 122/165 electrodes and higher grey matter-white matter (GM-WM) ratio for 99/165 electrodes. The computational efficiency of our algorithm enables near real-time planning of electrode implantations.
In this manuscript we focused on the development of our MTP algorithm leveraging existing methods for extracting the skull template, critical structures, and GM. Our algorithm was integrated into the EpiNav TM software platform to enable manual assessment of calculated trajectories. A larger, prospective, comprehensive clinical study of EpiNav TM is necessary to evaluate the utility of the software in planning intracerebral electrode implantations.