Undergraduate Research Assistant Brigham Young University, United States
Introduction:: Focused Ultrasound (FUS) is a non-invasive surgical technique that treats diseased tissues by thermal ablation. However, healthy tissues can suffer thermal damage during FUS surgery from heat transfer or stray sonication. During treatment planning, thermal modeling is used to predict tissue temperatures to avoid such damage. One method for thermal modeling discretizes the tissue in space and time and uses a forward-difference scheme to explicitly predict tissue temperatures during treatment. These explicit solvers suffer from numerical instability if the time steps are too large, requiring more computational resources to model treatments. A backward-difference scheme that implicitly predicts tissue temperatures would have no such limitation, enabling larger time steps and potentially reducing required computational resources. The purpose of this study was to modify an explicit thermal solver to be an implicit thermal solver and compare computation times for the same sonication and tissue conditions.
Materials and Methods:: A two-dimensional explicit solver of the Pennes bioheat transfer equation was modified to an implicit solver. Both solvers are MATLAB functions with inputs including a segmented tissue model, tissue properties, the sonication pattern over time, time step size, etc.
The implicit solver was verified by simulating a 60-second period with 30 seconds of constant power heating in the center of the model and 30 seconds of cooling. The tissue model was 10´10 cm2, segmented to 0.5-mm isotropic resolution with 3 randomly assigned tissue types with tissue properties comparable to those of human tissue (see Table 1).
To compare computation times, two test cases were simulated with inputs similar to the verification case. However, the sonication intensity was nonzero in every voxel with randomly assigned values. The first test case simulated 60 seconds with 30-second periods of constant sonication and cooling. The second test case simulated 3000 seconds with 1500 seconds of sonication and 1500 seconds of cooling. The sonication period included 50 random sonication patterns, each lasting for 30 seconds. For both test cases, the explicit solver was run at the maximum time step allowable for numerical stability, 0.25 seconds. The implicit solver was run with identical inputs. Next, the implicit solver was run with 30-second time steps, the maximum size that can model individual sonication periods. If the temperature results for the two implicit runs deviated more than 2% or 0.5 °C, the script was run with the time step halved, repeating until the results converged.
Results, Conclusions, and Discussions:: Results
The model verification simulation demonstrated that predicted temperatures agreed within machine precision for the explicit and implicit solvers.
The two test cases were run 10 times each to determine average computational runtimes. Table 2 summarizes the results for computation times. Figures 1 and 2 show the predicted temperatures for a single central voxel over time for both cases. For the first test case, the implicit solver was able to match the accuracy within 2% and 0.5 °C using 7.5-s time steps. For the second test case, the implicit solver was able to match the accuracy using 30-s time steps.
Conclusions
The implicit solver took significantly longer than the explicit solver for the same 0.25-s time steps. Simulating a 60-second period, the implicit solver was slower even with larger time steps. The majority of the implicit runtime is spent creating the coefficient matrix necessary for temperature calculations. As this step is performed only once, using the implicit solver to simulate longer time intervals alleviates this disadvantage. This is supported by the runtimes of the solvers in the second test case. There was a negligible difference in runtimes between the solvers when the implicit solver was run at larger time steps. Therefore, an implicit solver scheme is best suited to simulating longer times with large time steps. It is remarkable that the implicit solver can accurately predict the bioheat transfer for the entire sonication period in a single time step.
Discussions
The implicit solver’s largest computational cost was solving a large linear system of equations to predict new temperatures. For the given tissue model, the temperature coefficient matrix contains 1.6 million elements. However, this matrix is sparse and regular with only 5 nonzero diagonals. It is likely that algorithms exist to swiftly solve these spares matrices, which could greatly reduce runtimes for the implicit solver. Our next efforts will be finding or developing and then incorporating such algorithms.
