Simulation pipeline

Evaluating plans

Once the plan setup is complete, the plan can be evaluated by clicking the Evaluate button in the top toolbar. This exports a treatment planning file and uploads this to k-Dispatch. The plan is then queued for execution on a remote high-performance computing (HPC) server (see cloud computing). The plan evaluation happens asynchronously. This means it is possible to close k-Plan while a plan is being evaluated, or submit multiple plans for evaluation at the same time.

When the scheduled plan execution begins, the plan is evaluted by the SEM using the pipeline illustrated below. There are three main steps: aberration correction, forward planning, and thermal simulation. The aberration correction step is skipped if using a single element or annular array transducer, or if using a multi-element array and Run aberration correction simulation is unchecked in the Settings tab. The thermal simulation is skipped if Run thermal simulation is unchecked in the Settings tab.

graph LR
  A[planning<br/>file]-->B[aberration<br/>correction]-->C[forward<br/>planning]-->D[thermal<br/>simulation]-->E[results<br/>file]

The simulations are conducted using the k-Wave toolbox. The computing resources needed for the simulation (number of cores, amount of memory) are automatically selected based on the size of the simulation domain and the sonication driving frequencies.

Converting the primary planning image to material properties

For all stages in the simulation pipeline, k-Plan maps the primary planning image to material properties using the following pipeline:

graph LR
  A[CT Image]-->|resample<br/>and apply<br/>calibration|B[Density Image]
  B-->|segment|C[Skull Segmentation]
  B-->|segment|D[Soft-Tissue Segmentation]
  B-->|segment|E[Background Segmentation]
  C-->|mapping|F[Skull Properties]
  D-->|fixed values|G[Soft Tissue Properties]
  E-->|fixed values|H[Water Properties]

The resampling step resamples the primary planning image based on the simulation domain size and the specified number of grid points per wavelength. The image is then converted to density using the CT calibration.

The segmentation step uses an intensity based segmentation to divide the image into three regions: bone, soft tissue (e.g., skin and brain), and background (outside the head). The density thresholds used are \(\ge\) 1150 kg/m³ (bone) and \(\ge\) 900 kg/m³ (soft tissue). The segmentation masks are then optimised using morphological opening and closing operations.

Within the bone, the CT HU are converted directly to density and sound speed, while constant values are used for the acoustic attenuation and thermal properties (a constant value for the skull density is also used for the thermal simulations). Within the soft tissue and background, constant values are used for all properties. The properties for the background medium are based on the properties of water at 30^oC. The properties used in k-Plan are outlined in the table below.

Tissue	Density [kg/m³]	Sound Speed [m/s]	Absorption Coefficient [dB/(cm.MHz^y)]	Absorption Power Law (y)	Thermal Conductivity [W/(m.K)]	Specific Heat [J/(kg.K)]
Background	996	1509	2.17e-3	2	0.6	4200
Soft Tissue (Brain)	1045	1550	0.59	1.2	0.5	3600
Bone	Mapped from HU (1850 for thermal simulations)	1.33\(\times\)Density + 166.7	13.3	1	0.32	1300

Calcifications

The material properties are calculated using only the largest connected component of the skull segmentation. This means that calcifications in the brain may not be captured in the material property conversions.

Air

Currently, k-Plan does not include air in the simulations.

Aberration correction

Aberration correction simulations are used to calculate the phases for each element of a multi-element array such that the array focuses on the specified target position. The material properties used in the simulation are mapped from the primary planning image.

For each sonication, a simulation is performed where a virtual point source transmitting at the driving frequency is placed at the target position. The relative phase of the pressure signal at each transducer element is then recorded. The conjugated (reversed) phases are used for the forward planning simulations.

Aberration correction is not performed for single-element or annular array transducers.

Forward planning

Forward planning simulations are used to calculate the acoustic pressure field within the simulation domain. The material properties used in the simulation are mapped from the primary planning image.

For each sonication, the transducer model is created and positioned relative to the primary planning image. The transducer is driven using the specified sonication parameters. The relative amplitude and phase of the pressure field across the simulation domain are then recorded.

In all cases, the transducer is driven with a sinusoidally varying (continuous wave) signal. If custom element parameters (amplitude and phase) have been set for a sonication in k-Plan, these are directly used for the forward planning simulation. Otherwise, the amplitude and phase parameters for each element of the transducer are calculated as follows:

Single-element transducer: Phase is set to zero. Amplitude is set to give the specified target amplitude in free-field (water).
Annular array: Element phases are set to give the specified focal distance. Amplitude is set to give the specified target amplitude in free-field (water).
Multi-element array: If aberration correction is used, the phases are set to the output from the aberration correction simulation. If not, the phases are set based on geometric beam forming in free-field to focus at the specified target position. The element amplitudes are initially set to 1. If Re-scale target pressure is selected, the simulation output is re-scaled to give the specified target amplitude at the target position.

k-Plan uses the acoustic model in k-Wave for the aberration correction and forward planning simulations. This solves a set of coupled partial differential equations that account for linear wave propagation through a heterogeneous and absorbing medium. The numerical scheme uses a Fourier collocation spectral method to compute the spatial gradients, and a dispersion-corrected finite-difference scheme to integrate in time. Full details of the acoustic model can be found here.

To model the transducer, k-Plan uses the kWaveArray class in k-Wave, which allows the definition of transducer models that remain accurate regardless of the resolution of the computational grid. Full details of the transducer model can be found here.

Steady state

The acoustic simulations in k-Plan are conducted in the time domain. The pressure amplitude and phase within the simulation domain are extracted from the time-varying pressure signal recorded over the last acoustic period of the simulation.

At the beginning of the simulation, the acoustic pressure field will vary as the waves propagate across the domain and reflect from the skull bone. Once multiple reflections have had time to travel across the domain, the pressure field will reach a steady state. This means the extracted amplitude and phase don't depend on the precise length of time used for the simulation.

The length of time used for the acoustic simulations is controlled by the Aberration correction grid traversals and Forward planning grid traversals options in the advanced simulation settings. These set the simulation time to the specified number of grid traversals across the grid diagonal at the background sound speed.

The number of grid traversals it takes to reach steady state will vary depending on the skull and transducer geometries. The default setting of two grid traversals is generally a reasonable setting when using a simulation domain that covers the whole skull.

Thermal simulation

Thermal simulations are used to calculated the maximum temperature reached within the simulation domain, and optionally, the thermal dose. The material properties used in the simulation are mapped from the primary planning image.

k-Plan uses the thermal model in k-Wave for the thermal simulations. This solves the bio-heat equation assuming a heterogeneous distribution of material properties. The heating due the acoustic waves is calculated using \(Q = \frac{\alpha p^2}{c_0 \rho_0}\), where \(Q\) is the volume rate of heat deposition, \(p\) is the amplitude of the steady-state acoustic pressure field calculated by the forward planning simulations, \(c_0\) is the sound speed, and \(\rho_0\) is the mass density. The reduction in heating due to blood perfusion is ignored.

The evolution of the temperature field over time is computed using the specified pulse timing parameters. When the transducer is on, the heat source is set to \(Q\) within the simulation. When the transducer is off, the heat source is set to 0. If a cooling time is specified, after the pulse timing regime is complete, the thermal calculation is continued for the specified cooling time with the heat source set to 0.

Assumptions

The simulations make the following assumptions:

The wave propagation is linear. This is generally a good assumption for the driving frequencies and target intensities used for transcranial ultrasound stimulation (TUS).
Shear waves in the skull can be ignored. This is generally a reasonable assumption provided the transducer is positioned so the acoustic beam is close to normal incidence on the skull. However, if the transducer is positioned at oblique angles relative to the skull, the simulation will not be accurate.
The acoustic properties do not change with temperature. This is generally a good assumption for the small temperature changes seen in TUS.
The transducer can be modelled as an ideal radiator with optimised parameters. This is generally a reasonable assumption, provided that the idealised transducer definition has been validated against experimental measurements.
The transducer position is correct.
The acoustic and thermal properties can be accurately mapped from the primary planning image. This depends strongly on the quality of the primary planning image and the accuracy of the corresponding CT calibration file.
There is no blood perfusion in the brain. This is a conservative assumption.

If any of these assumptions are not met, then the simulations may not be correct.

Last update: October 4, 2024