Optimization of tumor treating fields using singular value decomposition and minimization of field anisotropy

Tumor treating fields (TTFields) are increasingly used to treat newly diagnosed and recurrent glioblastoma (GBM). Recently, the authors proposed a new and comprehensive method for efficacy estimation based on singular value decomposition of the sequential field distributions. The method accounts for all efficacy parameters known to affect anti-cancer efficacy of TTFields, i.e. intensity, exposure time, and spatial field correlation. In this paper, we describe a further development, which enables individual optimization of the TTFields activation cycle. The method calculates the optimal device settings to obtain a desired average field intensity in the tumor, while minimizing unwanted field correlation. Finite element (FE) methods were used to estimate the electrical field distribution in the head. The computational head model was based on MRI data from a GBM patient. Sequential field vectors were post-processed using singular value decomposition. A linear transformation was applied to the resulting field matrix to reduce fractional anisotropy (FA) of the principal field components in the tumor. Results were computed for four realistic transducer array layouts. The optimization method significantly reduced FA and maintained the average field intensity in the tumor. The algorithm produced linear gain factors to be applied to the transducer array pairs producing the sequential fields. FA minimization was associated with an increase in total current delivered through the head during a activation cycle. Minimized FA can be obtained for an unchanged total current level, albeit with a reduction in average field intensity. We present an algorithm for optimization of the TTFields activation cycle settings. The method can be used to minimize the spatial correlation between sequential TTFields, while adjusting the total current level and mean field intensity to a desired level. Future studies are needed to validate clinical impact and assess sensitivity towards model parameters.

each pair is activated in a 'square waveform' pattern. During activation, the individual sources deliver sinusoidal field waveforms at 200 kHz (for GBM). The induced fields are applied in approximately orthogonal directions based on a perpendicular layout of the two array pairs. The therapeutic efficacy of TTFields depends on the induced field intensity as well as the exposure time and the direction of the sequential fields relative to each other (Kirson et al 2004). Specifically, the growth rate of cancer cell cultures has been shown to decrease with increasing intensity and exposure-time of TTFields (Kirson et al 2004). Similarly, a recent study showed that the overall survival of GBM patients treated with TTFields correlated positively with the cumulative TTFields exposure, calculated as the average field intensity over time (Urman et al 2017). TTFields has also been shown to preferentially disrupt mitosis for cells dividing in the direction of the applied field (Kirson et al 2004, Wenger et al 2015a. This latter notion is supported by the in vivo results showing that two sequential and orthogonal fields increase the therapeutic efficacy by ≈20% compared to one constantly active field (Kirson et al 2007). This illustrates that the use of multiple field directions enables the distribution of the anti-cancer effect more evenly among cells dividing in different (random) directions thereby reducing the average tendency to differentially inhibit cells dividing in particular directions. Correspondingly, clinical TTFields therapy (Optune ® , Novocure, Ltd.) uses two pairs of transducer arrays activated in even 50% activation cycles, as explained above. The two corresponding fields cover two spatial directions and thereby target a larger fraction of cells compared to one field alone, although three linearly independent fields (ideally orthogonal and of equal magnitude) would be required to obtain an unbiased distribution. In a clinical treatment setting, the transducer array layout is planned to maximize the intensity of TTFields in the tumor (Trusheim et al 2017). However, in a recent study, we demonstrated that the complex anatomy and conductivity distribution of the head gives rise to considerable unwanted spatial correlation of the induced fields indexed as fractional anisotropy (FA). FA is extensively used in diffusion tensor imaging (Basser and Pierpaoli 2011). In the context of TTFields, it is calculated from the principal component magnitudes and measures the fractional deviation of the induced fields from the isotropic condition in which the sequential TTFields are perfectly orthogonal and have equal magnitude. FA generally takes values between zero and one, with higher values representing unwanted field correlation, i.e. the average field over an activation cycle will tend to have a directional preference of activity. Using the FA index, we have demonstrated considerable unwanted field correlation in the tumor region even for layouts planned to maximize treatment efficacy . This means that although the array layouts were planned to induce orthogonal (i.e. uncorrelated) sequential fields, the actual induced fields were in fact highly correlated (figure 1) to each other. This problem has previously been ignored in TTFields treatment and computational studies of TTFields potentially compromising or reducing the therapeutic efficacy.
In this study, we propose a computational method, which optimizes the activation cycle settings of TTFields to minimize unwanted spatial field correlation in any local region of a computational volume conductor model. The optimization method is based on a linear transformation and can be used for any particular array layout, any number of active fields, and any choice of desired average field intensity in the tumor. We provide modeling data validating the method's ability to decorrelate the induced fields for different layouts. The method may potentially be used in future embodiments of the TTFields technology to optimize the treatment efficacy for individual cancer patients.

Model preparation and electrical field calculations
The head model was created from MRI data obtained from a patient with left-sided parietal GBM (figure 2(A)). Written consent was obtained from the patient. The head model was segmented into six tissue types, namely, skin, skull, cerebrospinal fluid (CSF), gray matter (GM), white matter (WM) and tumor. The configuration of TTFields was equivalent to the Optune ® technology, which is used for clinical treatment. Four different array layouts were investigated (figure 2(C)). The distribution of the induced electric field was calculated using a finite element (FE) approximation of Laplace's equation for the electric potential, ∇ · (σ∇ϕ) = 0 (Miranda et al 2014, Korshoej et al 2016, 2017a, 2017b, Bomzon et al 2016, Lok et al 2017. This is valid at the low-to-intermediate frequencies of TTFields (e.g. 200 kHz), because the electromagnetic wavelength in the relevant tissues is much larger than the size of the head (Plonsey and Heppner 1967, Miranda et al 2014, Wenger et al 2015b. Furthermore, skin effect can be neglected at 200 kHz. The currents induced by TTFields are, therefore, mainly resistive (Ohmic) currents. Regarding the relationship between field frequency and cancer growth inhibition, we kindly refer the reader to Wenger et al (2015) and Kirson et al (2004). Finite element calculations were performed with SimNIBS (www.simnibs.org) (Windhoff et al 2013). We used Dirichlet boundary conditions with a fixed electric potential at each array (Saturnino et al 2015, Opitz et al 2015. Isotropic conductivity values were assumed for skin (0.25 S m −1 ), skull (0.010 S m −1 ) and CSF (1.625 S m −1 ) (Gabriel et al 2009, Korshoej et al 2016. For GM, WM and tumor tissues, we assumed anisotropic conductivity tensors inferred from diffusion MRI data (Tuch et al 2001, Korshoej et al 2016. The electric field was calculated as the numerical gradient of the potential. We calculated the current density using Ohm's law. The electric potentials, fields, and current densities were rescaled to obtain a peak-to-peak current amplitude of 1.8 A for each array pair. Further details about the head model generation and electric field calculation can be found in Korshoej et al (2016Korshoej et al ( , 2017a and Wenger et al (2018).

Estimation of average intensity and field correlation using singular value decomposition
For each finite element we defined the field matrixε with the transposed sequential field vectors E i in each row: (1) In addition, we defined the relative activation time α i of E i as where t i 0 is the 'on-time' of E i during the activation-cycle (α 1 = α 2 = 1 2 for Optune ® , i.e. a 50% activation cycle). In this context, it is important to appreciate, that the notion of activation cycle in this context refers to the pattern of field source activation and not the duty-cycle of the 200 kHz sinusoidal signal induced by the source. We used A R n×n as the diagonal activation time matrix with entries a ii = √ α i and define an 'activation-timeweighted' electric field matrix We performed a singular value decomposition (SVD) of P to obtain a new representation of P defined by up to three principal components, i.e. orthonormal (uncorrelated) basis vectors collected in the matrix W R 3×3 : The matrix Σ R n×3 contains the singular values σ k , k 3, i.e. the magnitudes of the principal components. We defined the average intensity of TTFields in each element as the Frobenius norm of P, i.e.
Thus, the average field was calculated as the square root of the activation-time weighted energy contributions from each field. The field correlation was estimated as the fractional anisotropy (FA) (Basser and Pierpaoli 2011), i.e.
FA estimates the fractional deviation of P from the isotropic condition with orthogonal fields of equal intensity.
In situations with <3 singular values the missing values were defined to be zero. For a more detailed discussion of the above derivations and definitions, the reader is kindly referred to Korshoej and Thielscher (2018). Figure 3 shows a schematic illustration of the calculation steps.

Optimizing the activation cycle of TTFields to minimize the fractional anisotropy
The principle component approach allows for optimization of the activation cycle such that FA can be minimized at a predefined average intensity for any array layout and in any region of interest. This can be done by applying the following linear transformation to P: Here s > 0 is the desired isotropic singular value, I nx3 R nx3 is the truncated identity matrix and Σ + R 3xn the pseudoinverse of Σ. In the following, we will refer to Q as the optimization matrix and we note that Q is symmetric. We see that the transformed field matrix G is given by which is clearly isotropic with two non-zero singular values equal to s as desired. Furthermore, G T = P T Q T = ε T AQ T , so we may think of the column vectors of G T as linear combinations of the original field vectors, in which the weighting factors defined by the entries in AQ T tell us how each field should be scaled to obtain an average distribution with minimum FA, i.e. equal singular values. To see this we define g j to be the j th row vector of G and thus g T j is the j th column vector of G T . Then, g T j = n i=1 q ji a ii E i and as expected we see that the scale factors q ji a ii of each field E i depend on both the corresponding activation time and the entries of the transformation matrix Q. Since the field vectors scale linearly with current, the j th row vector q j of Q may be considered a vector of relative gain factors q ji applied in the j th activation period of the electrode pair inducing the corresponding field. This gain may be considered a scale factor for a ii (activation time) or for the induced electric current, i.e. reduced anisotropy may be obtained by scaling current impressed by either layout or by scaling their relative activation times. To demonstrate how this linear transformation is able to significantly reduce unwanted FA of the time-weighted field distributions, we first need to define the region of interest (ROI) to which the optimization is applied. For a ROI comprising the finite elements k = 1, . . . , m we define the average optimization matrix Q as the median of optimization matrices Q k of all m elements in the ROI weighted by the relative volumes v k . If we consider the entries of Q as containing scale factors for the induced electric current density, then the field distribution of the optimized (on average) activation cycle is given by G =QP for each element, because the electric field depends linearly on the current density.

Results
In the following, we will provide examples of activation cycle optimization using four different array layouts. In addition, we will choose the desired isotropic singular value s so that E avr is maintained, i.e. we seek to derive scale factors for current (or activation cycle period) settings of the TTFields device, which maintain the same average intensity, albeit with minimized FA in the tumor region. The objective is to distribute the energy isotropically in the span of the induced field vectors (rather than across all three dimensions of physical space), since we are using only two array layouts in accordance with the standard TTFields application. This corresponds to the best achievable activation cycle optimization possible with the number of active transducer arrays (i.e. fields). To achieve this we will set s = , so that the average intensity E avr is maintained after activation cycle optimization. The potential impact of the optimization procedure is illustrated in figure 4, which shows a scatter plot of the singular values in each element of the field distribution in the tumor and the peritumoral regions before and after activation cycle optimization. The results were based on the field distribution induced by Layout 1 using standard current settings (1800 mA peak-to-peak) and relative activation times α 1 = α 2 = 1 2 . The figure confirms that the optimization procedure was able to reduce the difference between minimum and maximum singular values and hence reduce the FA in the ROI. Figure 5 shows the cumulative distribution functions of E avr (figures 5(A) and (C)) and FA (figures 5(B) and (D)) for four different layouts. Results are shown both before (solid) and after (stippled) activation cycle optimization. The optimization procedure maintains the E avr distributions largely unaffected, while FA is reduced for all layouts. FA was particularly high for Layout 1, which also produced the highest median E avr , and for this layout the optimization resulted in a pronounced reduction of FA. Thus the procedure was able to reduce FA and maintain mean intensity for all layouts, as expected, and the extent of FA reduction was higher for Layout 1, which had the most significant initial anisotropy.
The findings above are further illustrated in table 1, which gives the relative changes in the area under the curve (AUC) values for the distribution functions shown in figure 5 for all layouts and regions of interest. AUC was nearly unaffected (1.2%-2.9%) for E avr in both the tumor and the peritumoral region, while the corresponding values for FA were considerably reduced (1.7%-43%).
The distribution of scale factors contained in the optimization matrix Q k is shown in figure 6, using    absolute row sums of Q are > 1, indicating that the induced current density (or activation cycle duration) required to maintain E avr and minimize FA is higher than the default current settings. This can potentially be addressed by considering a normalized optimization matrix Q norm with entries q norm ij =q ij rank(ε) j=1 |qij| , i.e. each entry normalized by the corresponding absolute row sum. This would ensure a constant total current density during each period of the activation cycle. The resulting field intensity would be lower than E avr while FA would remain unchanged and minimized since the induced fields have the same direction but are scaled by the same factor. Consequently, the optimization process can be designed to consider the desired total current density, average field intensity and minimized FA in order to obtain an appropriate balance between these parameters.

Discussion
We have described a new method to optimize the activation cycle of TTFields for any given layout. The method enables minimization of unwanted spatial correlation (FA) of sequentially induced TTFields while maintaining the average field intensity at a desired level in a region of interest, or alternatively maintaining the total impressed current density at a particular level. The approach is based on a linear transformation of the distribution of principal field components, as derived by the authors in a separate study. The post-processing procedure is applied to calculated sequential field distributions in a volume conductor model. In this study, we have demonstrated that the algorithm robustly minimizes FA while maintaining the average field intensity for four different layouts. The optimization method produces a matrix of linear scale factors representing either (1) electrical current gains applied to each sequential field during the activation cycle, or alternatively (2) scaling factors used to balance the duration of each period in the activation cycle. If the optimization matrix entries are interpreted as current gain factors, then the proposed method uses a balanced simultaneous activation of all array layouts to produce resulting field vectors, which on average are less spatially correlated in any small volume in a region of interest such as the tumor. This is contrary to the contemporary implementation TTFields, which uses strictly sequential fields. When the entries of the optimization matrix are interpreted as modifying factors for the activation cycle duration, the fields would still be applied sequentially, although each period would be scaled by the corresponding entry. When the optimization matrix is used to scale the current, the reduced FA and unchanged field intensity would imply a higher total current to be impressed by the system, compared to the default settings, which might potentially increase the risk of tissue heating. The user would therefore be required to personalize the activation cycle settings to obtain a desired balance between total current, average field and FA in the tumor. When the optimization matrix is used to scale the activation cycle periods, the induced peak current density would be constant, however, a potential complication might be that FA minimization would imply longer relative activation time of the weaker field, while the stronger field would have a shorter relative activation time.
Despite the fact that the weaker field is active for longer periods, the difference in activation time may potentially affect the tendency to develop skin rash underneath the transducer arrays, which is one of the more commonly observed complications to TTFields treatment. Also, it is important to note that the proposed optimization approach may be affected by non-linearities in the relationship between the field intensity and the anti-tumor efficacy. Until now, the intensity/response relationship has only been properly characterized for the field-range between 60 and 240 V m −1 in vitro, in which the relationship is approximately linear (Kirson et al 2004). However, it is possible that there exists a lower threshold for anti-tumor activity (<60 V m −1 ), i.e. a field value below which no anti-tumor activity would be observed, which would imply that it would not be desirable to spend any time during the activation cycle at or below this level. Contrary, FA optimization would be more relevant in field ranges where the anti-tumor effect increases less than linear with increasing field strength. In those regimes, FA could potentially be reduced while largely maintaining the effectiveness of the array pair creating the peak field strengths. This could for instance be relevant if the anti-tumor activity would saturate at some higher field intensity (e.g. >240 V m −1 ). In conclusion, future studies are required to establish the dose/response relationship of TTFields more accurately and it is necessary to incorporate information about field intensity and correlation in such studies.

Limitations and future perspectives
A key ability of the proposed method is to reduce unwanted spatial field correlation, which is known to influence the extent of cancer cell damage in vitro. However, future studies are needed to validate the clinical importance of FA as it was recently done for field intensity (Urman et al 2017). In addition, it is necessary to characterize the independent correlation of each parameter with clinical outcome to be able to predict the expected treatment efficacy and produce accurate optimization algorithms. Also, it would be valuable to characterize the sensitivity of the FA and E avr estimates towards variations in model parameters (e.g. tissue conductivity) and tumor characteristics (e.g. morphology and location). Finally, it would be highly interesting to investigate potential performance limitations of the optimization algorithm, e.g. cases with very high FA in the region of interest. This might potentially be observed in cases with very high conductivity in the tumor or its vicinity due to shunting effects that would expectedly cause a locally inhomogeneous field distribution (Korshoej et al 2017a).