kWaveDiffusion

Time-domain simulation of heat diffusion and perfusion.

Syntax

kdiff = kWaveDiffusion(kgrid, medium, source)
kdiff = kWaveDiffusion(kgrid, medium, source, sensor)
kdiff = kWaveDiffusion(kgrid, medium, source, sensor, ...)
kdiff = kWaveDiffusion(kgrid, medium, source, [], ...)

Description

kWaveDiffusion is a class definition for the time-domain solution of the diffusion equation or Pennes' bioheat equation in 1D, 2D, and 3D. In addition to heat diffusion, Pennes' bioheat equation accounts for advective heat loss due to tissue perfusion (the flow of blood through tissue), and heat deposition (e.g., due to ultrasound absorption) [1]. The computation is based on a k-space pseudospectral scheme in which spatial gradients are calculated using the Fourier collocation spectral method, and temporal gradients are calculated using a k-space corrected finite difference scheme. For a homogeneous medium, the formulation is exact and unconditionally stable. For a heterogeneous medium, the time scheme allows larger time-steps to be taken for the same level of accuracy compared to conventional pseudospectral time-domain methods.

The most general partial differential equation solved by the kWaveDiffusion class is given by

A * dT/dt = div(Kt * grad(T)) - B * (T - Ta) + Q

where

A  = density [kg/m^3] * specific heat capacity [J/(kg.K)]
Kt = thermal conductivity [W/(m.K)]
B  = blood density [kg/m^3] * blood specific heat [J/(kg.K)] * blood perfusion rate [1/s]    
Ta = arterial temperature (blood ambient temperature) [degC]
Q  = volume rate of heat deposition [W/m^3]

In a homogeneous medium, the thermal coefficients are related to the diffusion coefficient (or thermal diffusivity) by

diffusion coefficient [m^2/s] = Kt / A

Note, when the diffusion coefficient is specified instead of the individual thermal coefficients, the equation that is solved is

dT/dt = div(D * grad(T))

For non-constant coefficients, this differs from the conventional heat equation (where the diffusion coefficient is taken outside the divergence operator). For convenience, the thermal coefficients related to perfusion can also be combined to give a single "perfusion coefficient" given by

perfusion coefficient [1/s] = B / A

The input parameters are assigned as fields to four input structures (kgrid, medium, source, and sensor) in the same way as the other models in the k-Wave toolbox. These structures define the properties of the computational grid, the distribution of medium properties, source terms, and the locations of the sensor points used to record the evolution of the temperature field over time.

The medium parameters can each be specified as a single scalar values in SI units (for homogeneous coefficients), or as matrices the same size as the computational grid (for heterogeneous coefficients).

The initial temperature distribution is specified by assigning a single scalar value or a matrix (the same size as the computational grid) to source.T0. A heat source can also be specified in the same way by defining source.Q (the volume rate of heat deposition).

The time history of the temperature field can be recorded automatically by specifying a series of sensor locations using sensor.mask. This is defined as a binary matrix (i.e., a matrix of 1's and 0's with the same dimensions as the computational grid) representing the grid points within the computational grid that will record the temperature field at each time step. The current sensor data can be queried at any point using the property kdiff.sensor_data (where kdiff is the name of the kWaveDiffusion object). The sensor_data is returned using MATLAB's standard column-wise linear matrix index ordering, indexed as sensor_data(sensor_point_index, time_index). If no time dependent output is required, the sensor input can be replaced with an empty array [].

After an object of the kWaveDiffusion class is created, the simulation is run by calling the method kdiff.takeTimeStep(Nt, dt), where kdiff is the object name, Nt is the number of time steps to take, and dt is the size of the time step. During the simulation, a visualisation of the temperature field is displayed. The current temperature can be queried (or modified) at any point using the property kdiff.T.

[1] Pennes, H. H. (1948). Analysis of tissue and arterial blood temperatures in the resting human forearm. Journal of Applied Physiology, 1(2), 93-122.

Inputs

kgrid	grid object returned by kWaveGrid
medium.diffusion_coeff	diffusion coefficient [m^2/s]
	OR
medium.density	tissue mass density [kg/m^3]
medium.specific_heat	tissue specific heat capacity [J/(kg.K)]
medium.thermal_conductivity	tissue thermal conductivity [W/(m.K)]
medium.perfusion_coeff	perfusion coefficient [1/s]
	OR
medium.blood_density	blood mass density [kg/m^3]
medium.blood_specific_heat	blood specific heat capacity [J/(kg.K)]
medium.blood_perfusion_rate	blood perfusion rate [1/s] (volumetric flow rate per unit volume)
medium.blood_ambient_temperature	ambient blood temperature within perfused tissue regions [degC]
medium.diffusion_coeff_ref	reference diffusion coefficient used within the k-space operator (default = 'max')
medium.perfusion_coeff_ref	reference perfusion coefficient used within the k-space operator (default = 'max')
source.T0	initial temperature distribution [degC]
source.Q	volume rate of heat deposition [W/m^3] (note, medium.density and medium.specific_heat must be defined when using source.Q)
sensor.mask	binary grid specifying where the temperature is recorded at each time step

Optional Inputs

Optional 'string', value pairs that may be used to modify the default computational settings.

Input	Valid Settings	Default	Description
'DisplayUpdates'	(Boolean scalar)	true	Boolean controlling whether details of the simulation are printed to the MATLAB command line.
'MovieArgs'	(string cell array)	{}	Settings for VideoWriter. Parameters must be given as {'param', value, ...} pairs within a cell array (default = {}), where 'param' corresponds to a writable property of a VideoWriter object.
'MovieName'	(string)	'date-time-kWaveDiffusion'	Name of the movie produced when 'RecordMovie' is set to true.
'MovieProfile'	(string)	'Uncompressed AVI'	Profile input passed to VideoWriter.
'PlotFreq'	(integer numeric scalar)	10	The number of iterations which must pass before the simulation plot is updated.
'PlotScale'	(numeric two element vector) or 'auto'	'auto'	[min, max] values used to control the scaling for imagesc (visualisation). If set to 'auto', a symmetric plot scale is chosen automatically for each plot frame.
'PlotSim'	(Boolean scalar)	true	Boolean controlling whether the simulation iterations are progressively plotted.
'RecordMovie'	(Boolean scalar)	false	Boolean controlling whether the displayed image frames are captured and stored as a movie using VideoWriter.

Outputs

kdiff

kWaveDiffusion object which can be used to run thermal simulations using the diffusion equation or Pennes bioheat equation

Dynamic Properties

Properties which can be queried or modified after the object is created.

Fieldname	Description
cem43	thermal dose given in cumulative equivalent minutes (cem) relative to T = 43 degC [mins]
T	current temperature field [degC]
Q	volume rate of heat deposition [W/m^3]

Static Properties

Properties which can be queried, but not modified, after the object is created.

Fieldname	Description
dt_limit	maximum time step for which the simulation remains stable [s]
lesion_map	binary matrix of cem43 >= 240 mins
lesion_size	total size of lesion_map (distance in 1D [m], area in 2D [m^2], volume in 3D [m^3])
sensor_data	time varying temperature recorded at the sensor positions given by sensor.mask [degC]

Methods

Methodname	Description
plotTemp	plot current temperature field in current figure window
setOptionalInputs('string', value, ...)	modify the optional inputs after the object is created
takeTimeStep(Nt, dt)	calculate the given number of time steps of the temperature field

See Also

bioheatExact