End-to-end workflow for finite element analysis of tumor treating fields in glioblastomas

Tumor Treating Fields (TTFields) therapy is an approved modality of treatment for glioblastoma. Patient anatomy-based finite element analysis (FEA) has the potential to reveal not only how these fields affect tumor control but also how to improve efficacy. While the automated tools for segmentation speed up the generation of FEA models, multi-step manual corrections are required, including removal of disconnected voxels, incorporation of unsegmented structures and the addition of 36 electrodes plus gel layers matching the TTFields transducers. Existing approaches are also not scalable for the high throughput analysis of large patient volumes. A semi-automated workflow was developed to prepare FEA models for TTFields mapping in the human brain. Magnetic resonance imaging (MRI) pre-processing, segmentation, electrode and gel placement, and post-processing were all automated. The material properties of each tissue were applied to their corresponding mask in silico using COMSOL Multiphysics (COMSOL, Burlington, MA, USA). The fidelity of the segmentations with and without post-processing was compared against the full semi-automated segmentation workflow approach using Dice coefficient analysis. The average relative differences for the electric fields generated by COMSOL were calculated in addition to observed differences in electric field-volume histograms. Furthermore, the mesh file formats in MPHTXT and NASTRAN were also compared using the differences in the electric field-volume histogram. The Dice coefficient was less for auto-segmentation without versus auto-segmentation with post-processing, indicating convergence on a manually corrected model. An existent but marginal relative difference of electric field maps from models with manual correction versus those without was identified, and a clear advantage of using the NASTRAN mesh file format was found. The software and workflow outlined in this article may be used to accelerate the investigation of TTFields in glioblastoma patients by facilitating the creation of FEA models derived from patient MRI datasets.


Introduction
Tumor Treating Fields (TTFields) therapy is an established treatment modality for glioblastoma, the deadliest form of malignant primary brain tumor (Pless et al 2013, Stupp et al 2015, Hottinger et al 2016. The alternating nature of the externally applied electric fields delivers penetrating energy at a frequency of 200 kHz through two pairs of transducer arrays placed onto the patient's scalp (Stupp et al 2015). Though several mechanisms of action have been identified, including disruption of the precise alignment and function of Tubulin and Septin during mitosis of tumor cells (Kirson et al 2004, 2007, Gera et al 2015, crucial facets of TTFields remain unknown. Specifically, how the strength of the electric fields, the surface structure of the tumor, and volume of cerebrospinal fluid affect treatment outcome remains unresolved. For example, there is no existing literature correlating electric field strength at the tumor with progression-free or overall survival. It is also unknown whether there is a 'dose threshold' that predicts better survival or what that field strength dosage may be. To answer these questions, patient progression-free survival and overall survival must be correlated against the electric fields experienced in the brains of patients, after accounting for confounding variables. Finite element analysis (FEA) on volumetric meshes representative of patient anatomy has the potential to reveal not only how TTFields are distributed but also how, if there is such a correlation between field strength and survival, to optimize their positioning to improved efficacy (Miranda et al 2014, Lok et al 2015, Wenger et al 2015, Huang et al 2017.
Creating FEA models for TTFields in the human brain generally involves tissue segmentation of anatomic data from magnetic resonance imaging (MRI) into 3D 'masks', the application of tissue dielectric properties for each finite element and computational analyses using finite element analysis software. The gold standard for volumetric mesh preparation has historically been through manual segmentation by an expert, but this comes with the obvious costs of time and effort, and the resulting segmentations are likely to be inconsistent between experts (Clarke et al 1995, Kaus et al 2001. Additionally, while a broad range of tools for segmenting images minimize some of these costs (Shen et al 2005, Mayer and Greenspan 2009, Ashburner 2012, Valverde et al 2015, creating FEA models for TTFields still requires time-consuming manual corrections, such as the removal of disconnected voxels, the incorporation of unassigned intracranial voxels, and the addition of 36 electrodes matching those of the TTFields transducers to the model's scalp. This approach is therefore not scalable to analyze a large number of patient datasets.
The difficulty in automated segmentation, cleanup, and electrode placement for FEA has prompted the development of two heuristic solutions. One involves scaling a pre-segmented patient brain to match the approximate dimensions of the patient MRI under analysis, with an ellipsoid representing the tumor (Hershkovich et al 2016), while another simplifies the patient anatomy into deformable convex hulls . But we have found in our previous works (Lok et al 2015(Lok et al , 2017 that small variations in patient anatomy, tumor shape, and cerebrospinal fluid volume can have profound implications for the predicted electric fields. We, therefore, developed a semi-automated workflow, using individualized patient anatomy and minimized user input, for the high-throughput computer modeling of TTFields. Our workflow consists of patient MRI pre-processing, co-registration, segmentation, and postprocessing using modified tissue probability maps. Included is the fully automated placement of electrodes, with conductive gels, on the surface of the model, mesh generation in ScanIP (Simpleware, Exeter, UK) and physics computation using COMSOL Multiphysics (COMSOL, Burlington, MA, USA). As expected, there was a greater degree of overlap between our workflow's end-product and the segmented masks with, as compared to those masks without, post-processing steps in ScanIP. Further, our workflow saves the considerable time required to generate personalized FEA models derived from patient MRI datasets and will be made available as a toolbox for SPM.

Methods
The procedure of our end-to-end workflow for FEA of TTFields is summarized in figure 1. Briefly, this process is divided into four distinct steps: (i) segmentation of MRI datasets using the New Segment algorithm from statistical parametric mapping 8 (SPM8) in MATLAB (MathWorks, Natick, MA), (ii) automated electrode and gel cylinder generation and placement, (iii) post-processing of segmented volumes using ScanIP and the addition of the gross tumor volume (GTV), manual correction and meshing, and (iv) physics simulation in COMSOL Multiphysics after defining all necessary parameters, material properties, boundary conditions and physics equations. MRI and pre-processing The MRI datasets were obtained based on an institutional review board-approved retrospective neuroimaging analysis of TTFields. Four patients, with ages ranging from 28 to 55 years, were selected based on differences in tumor geometry, size, location and response to treatment to ensure broad applicability of our workflow procedures. Response was classified by The Revised Assessment in Neuro-Oncology criteria (Wen et al 2010). The tumors of two responders were in the left frontal and right frontotemporal lobes (figures 2(a) and (b)), while the two non-responders had their tumors in the right frontal and right temporoparietal lobes (figures 2(c) and (d)). Only one of the patients did not have a necrotic core (figure 2(d)). For each dataset, a series of 156 contrast-enhanced sagittal slices, with a slice thickness of one mm, was acquired with 3D MP-RAGE (TE = 2.0 ms, TR = 8.4 ms) on a 3.0 Tesla GE HDxt scanner (General Electric, Fairfield, CT). The image matrix size was 512 × 512 pixels with a reconstructed resolution of 0.4883 mm by 0.4883 mm and a bit depth of 12. Furthermore, for contourer delineation of the titanium wires, GTV and necrotic core in ScanIP, T1 and T2 post-gadolinium image sets were captured with 24 images at a slice thickness of 5 mm at a resolution of 0.4883 mm by 0.4883 mm. All images were stored in standard Digital Imaging and Communications in Medicine (DICOM) format. For SPM8 implementation, electrode placement, and mesh generation in ScanIP, the computer used had two Intel Xeon E5-2630v3 CPUs (2.4 GHz) with 64 GB DDR3 RAM (2133 MHz).
The DICOM repository for patient image datasets is the sole input, where the images are first converted to Neuroimaging Informatics Technology Initiative (NIFTI) image format with SPM8's conversion function and then co-registered to a template in Montreal Neurological Institute space. The default workflow results in ten binarized masks: white matter, gray matter, cerebellum/brainstem, cerebrospinal fluid, orbits, skull, scalp, gel, electrodes and air. While several of these steps are common to the MRI segmentation workflow for various brain simulation studies (Huang et al 2013), the following sections highlight steps in which our workflow has been tailored specifically for the FEA of TTFields.
Tissue probability maps SPM8's New Segment algorithm, and the Unified Segmentation algorithm before it (Ashburner 2012), segment MRI datasets by estimating the likelihood of each voxel belonging to a particular tissue type based upon its intensity and position relative to other voxels assumed to be of the same tissue type (Ashburner and Friston 2005). Tissue probability maps encode spatial distributions for tissues a priori where the probability of a voxel in a given position belongs to a map's tissue is between zero (black) and one (white) after normalization of the map to the MRI (figure 3) (Ashburner and Friston 1997). In the Unified Segmentation algorithm, the concurrent segmentation and normalization of tissue probability maps are repeated until convergence of the voxel classification probabilities is reached (Ashburner and Friston 2005). The two non-responders had their tumors located in the right frontal (c) and right temporoparietal (d) lobes. Only one of the tumors did not have a necrotic core (d).
SPM8's New Segment algorithm segments by default the orbits as cerebrospinal fluid, but this may erroneously cause underestimation of electric field strength in the cerebrospinal fluid by inadvertently increasing the total volume of the fluid. To avoid this issue, we created another tissue probability map for the orbits and subtracted it from the tissue probability map of cerebrospinal fluid. Additionally, the extended tissue probability maps presented by Huang et al (2013) were expanded to incorporate the cerebellum and brainstem from Diedrichsen et al (2009).

Generation of electrode and gel masks
SPM's New Segment algorithm is an iterative model that cycles between registration of the patient images with the tissue probability maps and classification of each voxel to one of the tissue classes (Ashburner and Friston 2005). For the registration step, New Segment uses a set of cosine transform basis functions to adjust the tissue probability maps to the target image immediately after classification (Ashburner andFriston 1999, Klein et al 2009). A transformation, or deformation, is applied to the patient image, the tissues are classified, and the cycle repeats (Ashburner and Friston 2005). After a successful segmentation routine, the cumulative movement of the patient image to the tissue probability maps can be expressed as a triplet of scalar fields (for the X, Y, and Z directions) where each field describes the spatial deflection of voxels necessary to align with the templates. These scalar fields describing the non-linear transformation of the target image to template image, as calculated by unified segmentation, are deformation fields and their inverse estimates a non-linear transformation to normalize the template image referenced to the patient image (Litvak et al 2011, Tsuzuki andDan 2014). In this sense, an inverse deformation field is a map of the difference between the templates and the patient image (Chen et al 2008, Litvak et al 2011. To automate the placement of electrodes on the scalp, we embedded a series of cubes (intensity = 1.0, dimensions = 3 × 3 × 3 voxels) along the surface of the scalp tissue probability map as an 'electrode map'. Because the electrode map corresponds to the same space as the tissue probability map, the inverse deformation from New Segment generates electrode 'islands' along the scalp mask that the electrode generation routine extracts as seed points (Litvak et al 2011, Tsuzuki andDan 2014).
The MATLAB routine, however, requires another coordinate to adjust the resulting electrode seeds to the position, namely the center of the head. We added 6 head 'centering' points (intensity = 0.8, dimensions = 3 × 3 × 3 voxels) to the 'electrode map' near the locations of the template's vertex, nasion, inion, right and left preauricular points, and the foramen magnum. The routine finds them in the deformed electrode map by their islands' maximum intensity; islands with maximum intensities greater than 0.95 are classified as electrode seed coordinates while islands with maximum intensities between 0.70 and 0.81 are classified as centering coordinates. The routine uses the mean of these 6 centering coordinates as the head's center.
Several sequential processes generate gel and electrode cylinders on the surface of the scalp. For each, the intensity is referenced to coordinate positions in the binarized scalp mask.
I. Electrode seed points are moved to the surface of the scalp (figure 4(a)). First, a center-toseed vector between the head's center to the electrode seed is created. If the intensity at the electrode seed's initial coordinate is equal to 'one', the seed point is repeatedly moved outward along the normal of the tangent to the surface of the scalp until the intensity at its next position along the center-to-seed vector would be 'zero'. If the seed's initial intensity is 'zero', an array of coordinates between the seed and the center is created. If none of the coordinates have an image intensity of 'one', the routine classifies the seed point as inside the scalp and vice versa. If the seed is inside the scalp, it is moved outwards until its intensity becomes 'one', and outwards further until the voxel before it reaches empty space. Seeds initially outside the scalp are moved inwards until they reach a coordinate with an intensity of 'one'. II. Electrode seed neighbors are adjusted to be approximately 22.5 mm apart from one another (figure 4(b)). A nearest-neighbor search finds surface seed neighbors. Seeds with two neighboring seeds less than 30 mm in distance are assumed to be the center of a triplet of electrodes, and the two neighboring seeds are moved towards or away from the center seed, along the vector between the two, if their distance from the central electrode seed is less than 21.5 mm or greater than 23.5 mm. The routine repeats Step I once more to move adjusted electrodes back to the surface of the scalp. III. The normal vector for the gel and electrode is determined by principal component analysis (figure 4(c)). A list of coordinates within a 10 mm radius sphere around the seed point is collected, and principal component analysis is carried out to find the normal vector to the surface of the scalp (Koessler et al 2008). This radius is twice the magnitude of the electrode's to overcome anomalies in the local surface structure. IV. The starting point of the gel, the intersection point between the gel and the electrode, and the end-point of the electrode are then specified (figure 4(d)). A circle of points in 3D space with a 10 mm diameter, perpendicular to the surface normal vector, is created at the seed point. The circle is moved towards the scalp's center until all points in the circle have an intensity of 'one', and its center is stored as the base point for the gel cylinder. The circle is then moved outwards along the normal vector until the first position at which all points in the circle have an intensity of 'zero'. Because the electrodes emitting TTFields sit on a layer of gel and not directly on the scalp, the circle is moved out an additional 2 mm to a point marked as the intersection of the electrode and the gel layer. Finally, the end point of the electrode cylinder is found 2 mm further outward, matching the height of the TTFields electrode.
With the start and end points calculated for both gel and electrodes, the routine generates gel and electrode cylinders with a diameter of 10 mm and saves them as binary masks in two NIFTI files. On average, the time required for tissue segmentation in our workflow requires 25 min for each patient image dataset and the electrode-generating routine requires 5 min.

Post-processing in ScanIP
By default, binary masks are saved in Analyze image data format for importation into ScanIP as background images, after which ScanIP masks are created by threshold filling. In FEA, well-defined boundary conditions are a major prerequisite for producing correct numer ical models. Additionally, the presence of island cavities, artifacts resulting from the surface architecture of a manually segmented tumor or other tissue masks after auto-segmentation, can make meshing problematic and/or unintentionally increase the number of degrees of freedom to be solved. To mitigate these potential pitfalls, several post-processing procedures are performed in ScanIP that clean the masks: I. A Gaussian filter function smooths the boundaries between masks to avoid convergence issues. II. Cavities in the mask are filled, and islands above a threshold (varies with tissue type) are removed. III. The current mask is duplicated and then dilated (by one to three voxels, depending on the tissue mask) and Boolean added to the 'next mask' on all slices.
To avoid issues of overlap, the post-processing script wraps tissue masks around one another in a concentric manner such that (i) gray matter wraps the dilated duplicate of white matter, (ii) cerebrospinal fluid wraps gray matter and cerebellum, (iii) the skull wraps cerebrospinal fluid, and (iv) the scalp wraps the skull and orbits (figure 5). This scripted post-processing requires approximately 6 min per image dataset.

Contour of tumor and meshing
T1 and T2 DICOM image sets are converted to NIFTI format in SPM8 and are registered to the NIFTI image that was input to New Segment. These registered images are saved in Analyze format and imported to ScanIP as background images to assist with segmentation. The attending physician segments the GTV as well as the necrotic core, and titanium wires if present, in ScanIP using built-in segmentation tools. Automated tools for this exist (Clark et al 1998, Bauer et al 2011, Havaei et al 2017, and more are under development (Menze et al 2015), but we have accepted the trade-off in time for the ensured accuracy of the FEA models from slight variations in patient anatomies. Time needed to segment the GTV and necrotic core can vary based on their respective sizes, but 1 h is common. However, an interpolation tool and threshold painting tool in ScanIP can be used to speed up the contouring process. Meshes were generated in ScanIP with the +FE Free mesh algorithm using linear tetrahedral geometry for elements and both MPHTXT and NASTRAN mesh formats as output (Yagawa 2011). Twenty iterations for smart mask smoothing were used during the mesh pre-processing procedure. A fully meshed model ready for export to COMSOL took approximately 30 min.

Finite element analysis physics preparation and computation
The finalized mesh file was first imported into COMSOL Multiphysics using the mesh import interface. Minimal boundary partitioning was selected as an option for both MPHTXT and NASTRAN file import. However, the MPHTXT file format requires the manual selection of tissue domains representing each of the masks, while the NASTRAN file format automatically imports pre-selected mask domains from the mesh data for the subsequent aggregate assignment of material properties in COMSOL. Electric conductivity and relative permittivity, reported in table 1, were calculated and assigned to white matter, gray matter, cerebrospinal fluid, skull, scalp, orbits and cerebellum/brainstem (table 1) after calculation for the properties at 200 kHz using the Gabriel dispersion relationships (Gabriel 1996, Hasgall et al 2015. Properties for GTV and necrotic core were found in the work of Peloso et al (1984); gel properties were calculated based on impedance values measured by the manufacturer of the gel used in combination with TTFields; electrode properties were found in Arlt et al (1985) and Wenger et al (2015). Titanium relative permittivity and conductivity were found in Wypych et al (2014) and in reported values by Eddy Current Technology Incorporated (Virginia Beach, VA, USA), respectively. The boundaries for each electrode in each of the four sets of arrays were explicitly assigned with a waveform function representing the sinusoidal input of the TTFields. Joule-Heating physics was integrated into each model to account for any effects from the heating of tissues through extended exposure to TTFields. Either frequency-dependent or transient study steps can be used to parametrically compute for the electric field distribution throughout the brain at discrete time points as explained by Lok and Sajo (2016). The time point associated with the root mean square (RMS) value was chosen for cross-model comparisons since the RMS of a sine wave is the arithmetic mean and appears four times per cycle, and use of the RMS time point provides a means to determine the average exposure to the TTFields over an extended period among different patients. The models used in the present study were computed on a custom-built computer consisting of an i7-4690X (4.1 GHz) with 64 GB of DDR3 RAM (2133 MHz). Each frequency-dependent model took approximately 30 min to compute using approximately 20-30 GB of RAM, while each transient model took approximately 5.5 h to compute using approximately 30-50 GB of RAM (both depending on the number of element domains, boundaries, edges, and vertices). An additional 12 h of derived value volume integration was required to evaluate the data necessary to construct electric field-volume histograms, including greater than 4 k, 600 DPI high resolution stacked image export of electric field distribution with intensity represented by either colors or grayscale (in axial, sagittal and coronal anatomical orientations) and with slice thickness varying between 1 and 3 mm, with or without electric field vector overlay.

Evaluation of workflow
To investigate the accuracy and precision of the electrode generation routine, we measured the distances between electrodes along the central long axis and three short axes of each electrode transducer. We found the center of each electrode in ScanIP and calculated the two Euclidean distances from the middle electrode to its two outer electrodes in each of the four specified axes. This process consists of 32 measurements for a single brain, eight per transducer, and 128 measurements for all four brains.
We also compared the results of segmentation-the binarized output from SPM8 New Segment-between models with post-processing versus models without post-processing, and models with post-processing and manual correction (as the baseline). Since an individual who performs the segmentations is typically familiar with neuroanatomy, an assessment was made between the manually corrected, semi-automated models versus the automated models with and without post-processing on a mask-by-mask basis. Because the manual segmentation of the GTV and necrotic core are integral to the workflow, we kept these two masks constant across models with respect to the patient, copying these two masks from the manual model to the other models.
The Dice coefficient (DC), or Sørensen index, is a common metric for measuring the degree of geometric similarity between two structures, with higher scores representing a greater degree of overlap. The Dice coefficient was used in this study to quantify the percentage overlap between the models for each patient on a per-tissue mask basis (Babalola et al 2009).
We calculated the average relative difference in electric field maps derived from automated models with post-processing versus semi-automated models using the following formula below where E AUTO is electric field strength at a given voxel from the automated models with post-processing while E SEMIAUTO is the electric field strength at the same voxel from a semiautomated model ( figure 7). As stated previously, the GTV as well as the necrotic core and titanium wires masks, if present, were copied across models for each patient for comparison sake. Due to discontinuous boundary conditions between masks in the automated models without post-processing, from unclassified space, electric field maps were not computed for these models. Reported relative differences of the fields are averages over all voxels within the region of the head in the electric field maps, excluding empty voxels outside the head.
Finally, to compare mesh file formats generated by MPHTXT and NASTRAN, the relative difference in electric field maps was calculated in the same manner described above, with NASTRAN as the reference map.

Results
Our workflow encompasses MRI pre-processing, co-registration, segmentation and postprocessing. After fully automated placement of electrodes and conductive gels on the surface of the model, a mesh is generated using ScanIP and results are computed using COMSOL Multiphysics. Automated segmentation is implemented with the SPM8 New Segment algorithm along with modified extended tissue probability maps to include cerebellum and orbits in addition to standard tissues such as white matter, gray matter, cerebrospinal fluid, skull and scalp. Further automated post-processing in ScanIP merges disconnected voxels, smooths surfaces, and layers the masks to prevent erroneous boundary conditions. Additionally, postprocessing assigned all unclassified intracranial voxels to the neighboring tissue types, most often cerebrospinal fluid (figure 6(c)). Despite the existence of tools for the automatic identification and segmentation of tumors from MRI, we opted for manual segmentation of both the GTV and the necrotic core, given the need for accuracy and precision in the delineation of heterogeneous tumor geometry, size, and location.
We investigated the results of this workflow in several ways. First, we evaluated the accuracy of our automated electrode generation method by comparing the distance of neighboring electrodes derived from our method to the actual measured distance on the transducer arrays. The average distance between neighboring electrodes along the short axis of each array was 23.5 ± 0.8 mm (range 21.8-25.5 mm) and along the central long axis was 45.6 ± 2.2 mm (range 42.0-50.3 mm). This is similar to the distance of 22.5 mm on the short axis and 45 mm on the central long axis as measured on the array, and the variance of our method is about 4% and 5%, respectively.
Second, using the four patient datasets, we applied the Dice coefficient by comparing the extent of spatial overlap across segmented masks from SPM8 New Segment without postprocessing in ScanIP (figure 6(b)), SPM8 New Segment with post-processing in ScanIP (automated model, figure 6(c)), and our full end-to-end workflow model (semi-automated model that incorporated manual correction, figure 6(d)). Comparing the segmented masks without and with post-processing to the full workflow model, the spatial overlaps between the six largest maps (scalp, skull, cerebrospinal fluid, gray matter, white matter, and cerebellum/ brainstem) increased from an average Dice coefficient of 70.7 ± 16.8 to 85.2 ± 10.0 (P = 0.005, twotailed paired student t-test on a mask-by-mask basis, figure 6). The tissue masks that exhibited the greatest improvement are cerebrospinal fluid and skull, with their Dice coefficients improved from 38.3 ± 5.2 to 65.3 ± 11.2 and 60.5 ± 10.9 to 81.4 ± 7.0, respectively. After post-processing in the automated model and compared to our semi-automated model, the mask for cerebrospinal fluid had the least overlap while the scalp had the most overlap, and their Dice coefficients were 65.3 ± 11.2 and 96.2 ± 1.4, respectively. Collectively, ScanIP post-processing of the segmented masks generated by SPM8 New Segment helps to construct the final mesh needed for the FEA of TTFields.
Third, we compared COMSOL generated electric field maps for automated models with post-processing versus semi-automated models. We found an average relative electric field difference of 5% ± 3.8% on a voxel-to-voxel basis between the automated with post-processing and semi-automated models (figure 7), suggesting a discernible difference as a result of the manual correction to the automatically generated models. In addition, the cumulative electric field intensity for each tissue, which is represented by the area under the curve of the electric field-volume histogram, was compared between the automated and semi-automated models (figures 7(d)-(f)) and an average relative difference of 12.2% ± 5.7% was found in all major masks.
Additionally, mesh file formats MPHTXT and NASTRAN were compared using a similar electric field relative difference analysis to investigate the potential difference in electric field maps because of mesh format choice. There was little to no significant difference and the averaged relative difference in electric field strength was 0.9% ± 1.0% by the comparison of the full maps ( figure 7(c)). Similarly, by relative difference comparison of electric fieldvolume histograms (figures 7(d) and (e)), little difference was noted with 3.5% ± 6.0% for scalp, 2.0% ± 3.4% for skull, 0.5% ± 0.6% for cerebrospinal fluid, 0.04% ± 0.04% for gray matter, 0.07% ± 0.06% for white matter, 0.04% ± 0.02% for GTV, 0.3% ± 0.3% for necrotic core, and 4.5% ± 7.5% for combined cerebellum/brainstem. Because the differences in FEA from using MPHTXT or NASTRAN file formats are minute, we have since incorporated NASTRAN as the standard mesh file format in our workflow. This is because the NASTRAN file format obviates the cumbersome and potentially error-prone procedure of manual selection of each domain and the assignment of material properties that are inherent in using the MPHTXT file format.
Finally, we noted a small difference in the average electric field strength between responders and non-responders, with the two responders exhibiting an average electric field strength of 68 V m −1 ± 17 V m −1 versus 58 V m −1 ± 5 V m −1 in non-responders. Further, the models of responders had higher average necrotic core electric field strengths at 41 ± 7 V m −1 versus 34 V m −1 in the one non-responder with a necrotic core. However, the sample size was too small for any statistically valid comparison.

Discussion
In this study, we presented an efficient and comprehensive semi-automated, end-to-end workflow for the high-throughput generation of personalized patient finite element models for the study of TTFields. This workflow incorporated modified, extended tissue probability maps along with an electrode generation routine written in MATLAB that uses inverse deformation fields from segmentation to automate the placement of electrodes on the scalp imitating the positions of the transducer arrays emitting TTFields. Furthermore, we present a method for post-processing FEA models in ScanIP and physics simulation using COMSOL Multiphysics that included model evaluation as an example of the analyses expedited by our procedures.
We processed MRI datasets from four patients using our workflow to investigate its quality and speed. The output of the automated portion of the workflow consists of binarized masks for white matter, gray matter, cerebrospinal fluid, brainstem/cerebellum, skull, orbits and scalp, as well as extra-cranial electrodes and gels. Compared to the final model used in our workflow, there was a significant improvement in the Dice coefficient of the segmented masks after post-processing by ScanIP, suggesting that this post-processing procedure is important for constructing the final mesh product needed for the FEA of TTFields. A physician familiar with patient neuroanatomy then adds the GTV and, if present, the necrotic core and titanium wires. Segmentation, electrode generation, post-processing, and mesh model generation took on average 25, 5, 6, and 30 min, respectively, and the entire automated portion of the workflow for mesh preparation took, in aggregate, slightly over an hour. But this timeframe may vary depending on the size and resolution of the MRI scans. Manual segmentation of the GTV, necrotic core, titanium wires and other minor adjustments of masks requires a variable length of time, depending on the complexity of these structures, but 1.5 h has been our experience. Total time for one fully computed and completely presentable personalized FEA model took approximately 24 h, not including idle time between completed steps. Therefore, our workflow can generate personalized FEA models derived from patient MRI for analysis of large patient datasets and, eventually, treatment planning in the real-world clinical setting.
We created a method for the automated placement of electrodes because the manual placement was cumbersome and prone to error (Clarke et al 1995). With manual electrode placement, one would need to add 72 cylinders in the correct locations, which consist of a total of 36 electrodes and their corresponding gel layers. We do not know of any other existing, expedient, and consistent solution for positioning cylindrical electrodes with a gel layer on the scalp. The closest ones are several semi-automated solutions for placing electroencephalogram Figure 7. Single 2D axial slice representation of electric field intensity. The electric field map from fully automated segmentation (a), the field map from semi-automatic segmentation with manual correction (b), and the relative difference between the two on a pixel-by-pixel basis for each slice (c). For aggregate comparisons, 3D reconstruction from all axial slices of the electric field distribution maps was co-registered, and the relative difference between the electric field distribution was computed on a voxel-byvoxel basis. Representative electric field-volume histograms for GTV (d) and necrotic core (e) compared between automated models with post-processing versus semiautomated models with manual correction. The averages for the relative difference for the electric field-volume histograms of each mask type were tabulated (f). electrodes (Oostenveld et al 2011, Huang et al 2013, but none position the electrodes flush against the scalp in pre-specified locations. The method used in our workflow results in an acceptable variance of 4% and 5% with respect to electrode placement along the short and central long axes of the transducer array, respectively. Still, our approach has two possible sources of error. First, there is a small inherit error factor associated with the transformation defined in the inverse deformation field (Chen et al 2008) and, second, the resulting distances between the electrodes along the surface of the scalp may vary across heads of significantly different sizes. However, the routine partially corrects for this error along the transducers' short axis, via a nearest-neighbor search with fine adjustment ( figure 3(b)), and we have not encountered a segmented mask that required significant adjustment.
To investigate segmentation results before and after scripted post-processing in ScanIP, we computed the Dice coefficient to evaluate the extent of the spatial overlap between the 6 largest masks (scalp, skull, cerebrospinal fluid, gray matter, white matter, and cerebellum/ brainstem). The average Dice coefficient increased significantly by approximately 21%, or from 70.7 to 85.2 ± 10.0, on a mask-by-mask basis. The mask for cerebrospinal fluid exhibited the greatest improvement after post-processing in ScanIP but it had the lowest absolute Dice coefficient or spatial overlap with the final workflow model. One possible explanation for this variance may be due to the manual addition of patients' tumors and titanium wires within the skull, both of which could affect the segmentation results (Crinion et al 2007, Klein et al 2009. Extended workflows with increased image resolution, reduction of MRI bias, and multi-modality image input could potentially lead to enhanced results. The average relative difference in electric field strength was 5% between the automated and our semi-automated workflow models, suggesting the existence of small variations between models after manual adjustments as compared to those that were fully automated. Furthermore, there was an average of 12.2% in relative difference for the electric field-volume histograms of the masks between automated models after post-processing compared to the semi-automated models. These data suggest that our workflow to generate electric field maps may differ in its results from manual adjustments when compared to the automatically generated models. Unfortunately, the extent of variation caused by manual correction is almost entirely dependent upon the contouring decisions made by the experts while they are adjusting the final workflow models. For example, the boundary between gray and white matter is rarely clear cut, and studies suggest that there is about 20% variance in the delineation of gliomas by the same expert, even among those who are experienced (Mazzara et al 2004). As such, it is difficult to objectively quantify the accuracy of the described workflow, even with manual correction, since the individual might display a degree of variance in their assignment of tissue type into the masks. Our approach in the workflow initially used meshes formatted in MPHTXT files, but these files created a number of uncertainties and assumptions since one had to manually select each domain, which was often complicated by domains that were split into a large number of smaller elemental domains. Manual selection also created questions of whether one was associating the right material properties with the correct mask or tissue type and this process was extremely cumbersome and time-consuming. We have since used NASTRAN as the standard mesh file format of choice. This file format automatically includes pre-selected domains directly from the models' masks prepared in ScanIP, and we have found little difference in the average electric field maps of models generated from either MPHTXT or NASTRAN file format.
A more relevant metric for quantifying the workflow's performance would be the degree to which its predicted electric field map correlates with those measured in the brain of a patient undergoing TTFields therapy. A related investigation was recently carried out for the analysis of transcranial electric stimulation in 10 patients and it was found that current flow models, generated by a finite element workflow not dissimilar from the one described here, accurately predicted those electric fields that were measured in vivo (Huang et al 2017). A similar intracranial investigation of TTFields might shed light on the predictive ability of FEM-based TTFields models.
We have noted a non-significant difference in the average electric field strength between responders and non-responders, with the responders exhibiting an average electric field strength of 68 ± 17 V m −1 versus 58 ± 5 V m −1 in non-responders. The models of responders had average necrotic core electric field strengths of 41 ± 7 V m −1 versus 34 V m −1 in the one non-responder with a necrotic core. This comparison is reported for demonstration purposes only, highlighting the potential output of the described workflow, but the number of patients is too low for such a survival analysis. To account for confounders such as age, the extent of surgical resection and electric field strength, a minimum of 15 patient models would be necessary (Vittinghoff and McCulloch 2007). Our workflow is designed for the high-throughput analysis of TTFields and is the first to have been characterized against patient MRI with real tumors, but other software has been presented to address similar goals. The first and least related is the NovoTAL ™ System, a treatment planning software that generates images of a transducer layout believed to administer optimally TTFields at the site of the glioblastoma (Chaudhry et al 2015) after input of 20 anatomical measurements by a certified physician. The software is not presently designed to produce models for the FEA of TTFields. The second proposed solution is the use of convex hulls, in which hulls approximating anatomical structures could be deformed to the anatomy of the patient after measurement of the patient scalp or, in addition, measurement of internal structures . This approach obviates the need for electrode placement, and it is easier to implement computationally. The third proposed solution, which is principally related to the one described previously, involves scaling a pre-segmented human head, an atlas, to match the patient anatomy under analysis (Hershkovich et al 2016). Strengths of this approach are the automated analysis and minimal software dependencies, such that the barrier to entry is low. But the idealization of the tumor as an ellipsoid, which is an approximation, would have failed for two of the four patient anatomies used here as test cases. Both patients had diffuse and non-spherical tumors that would not have been accurately represented in the prototype solution proposed in Hershkovich et al (2016). Furthermore, stretching a pre-segmented atlas brain based on manual measurement from the user is imprecise, particularly given the existing registration algorithms that automate the task (Klein et al 2009). Finally, as we have found that small changes in the cerebrospinal fluid volume altered the strength of TTFields delivered to the GTV (Lok et al 2017), both the second and third software solution for FEA of TTFields are liable to produce field maps that vary from those that result for patient-specific MRI segmentation. Unified segmentation, which is inherently more granular to the patient anatomy being investigated (Ashburner andFriston 2005, Klein et al 2009), can be utilized for FEA with ease after accounting for auxiliary concerns such as pre-processing, post-processing, and electrode placement, as described in our work, and we believe this approach gives a more accurate view of TTFields within the intracranial space of patients.

Conclusion
Using real and individual patient MRI datasets, we described the first end-to-end, semiautomatic segmentation-based workflow procedure for generating personalized finite element models for the delineation of intracranial TTFields. MRI pre-processing, segmentation, postprocessing, electrode placement and physics computation are described such that the procedure is mostly automated and can be applied in a timely fashion for use in real-world clinical settings. Only the addition of GTV, necrotic core and Titanium wires requires manual input. We believe that this workflow will facilitate inquiries into the fundamental physics of how TTFields propagate in the human head and how they might be optimized to improve therapeutic benefit against glioblastoma.