Theoretical Modeling of the Bioheat Equation Using Greens Theorem


Bioheat transfer processes in living tissues are often influenced by the influence of blood perfusion through the vascular network on the local temperature distribution. When there is a significant difference between the temperature of the blood and the tissue through which it flows, convective heat transport will occur, altering the temperatures of both the blood and the tissue. Perfusion based heat transfer interaction is critical to a number of physiological processes such as thermoregulation and inflammation. The blood/tissue thermal interaction is a function of several parameters including the rate of perfusion and the vascular anatomy, which vary widely among the different tissues, organs of the body, and pathology. The literature contains an extensive compilation of perfusion rate data for many tissues and organs. The rate of perfusion of blood through different tissues and organs varies depending on factors such as physical activity, physiological stimulus and environmental conditions. Further, many disease processes are characterized by alterations in blood perfusion, and some therapeutic interventions result in either an increase or decrease in blood flow in a target tissue.

Bagaria, and Johnson (Bagaria, and Johnson, 2003, 2005) studied the the bioheat equation numerically and analytically for hyperthermia applications of cancer treatments. The model studies the proper distribution of magnetic particles throughout the tumor could minimize the damage to the surrounding healthy tissue while still maintaining a therapeutic temperature in the tumor. However, this distribution is defined mathematically but is not feasible to control in practical applications. And the analytical solution is complicate and obscure. Rosensweig, (Rosensweig, 2002) present a model for heating magnetic fluid with alternating magnetic field.

In this chapter, we the bioheat equation with a heat source of constant density power located as the center of a spherical domain using Greens functions. The point heat source model the heat generated by magnetic particles under the effect of an alternating magnetic field, used in studies of local magnetic hyperthermia. Parametric studies of the temperature profiles are carried out to study the effect of different parameters like the heat generation rate, perfusion rate and diameter of the point source on the maximum temperature and on the temperature profile. Some discussion about important parameters research issues in cancer hypertermia are addressed.

Mathematical Formulation

Pennes model (Pennes, 1948) describes the effects of metabolism and blood perfusion on the energy balance within tissue. Basically, these two effects are incorporated into the standard thermal diffusion equation and the equation is called the bioheat equation. Here, the domain of study will be a spherical domain with a heat source of radius \(r_1\) at the center of the spherical domain of radius \(R_0\) as is shown schematically in the figure below.

Model for an internal heat source inside of a tissue.

The bioheat equation in spherical coordinates with an internal heat generation at the center of the sphere can be written as:

\(\begin{align} \large \frac{\rho C_p}{k}\frac{\partial T}{\partial t} = \frac{1}{r^2} \frac{\partial}{\partial r}\Bigg(r^2\frac{\partial T}{\partial r}\Bigg) + \dot{q_p} + \frac{\dot{q_m}}{k} + \frac{\dot{Q(r)}}{k} \end{align} \)

where \(\rho\), \(C_p\), \(k\) are the density, specific heat and thermal conductivity of the tissue, \(\rho_b\), \(C_{pb}\) are the density and specific heat of the blood, and \(q_m\) is the metabolic heat generation. Penne’s model for the perfusion heat source is:

\(\begin{align} \large \dot{q_p} = \frac{\omega\rho_b C_{pb}}{k}(T_{\infty} - T) \end{align} \),

where \(\omega\) is the perfusion rate (\(m^3/s\) of volumetric blood flow per \(m^3\) of tissue), \(T_{\infty}\) is the arterial temperature and \(T\) is the local tissue temperature. Then, by incorporating the Penne’s model into the diffusion equation (1) yields

\(\begin{align} \large \frac{\rho C_p}{k}\frac{\partial T}{\partial t} = \frac{1}{r^2} \frac{\partial}{\partial r}\Bigg(r^2\frac{\partial T}{\partial r}\Bigg) + \frac{\omega\rho_b C_{pb}}{k}(T_{\infty} - T) + \frac{\dot{q_m}}{k} + \frac{\dot{Q(r)}}{k} \end{align} \)

We will define the following dimensionless parameters:

\(\begin{align} \large r^* = \frac{r}{R_0}, t^* = \frac{\alpha t}{R_0^2}, T^* = \frac{T - T_{\infty}}{T_{\infty}} \end{align} \)

Writing equation 3 using these dimensionless parameters, we get the following:

\(\begin{align} \large \frac{1}{r^{*2}}\frac{\partial}{\partial r^{*}}\Bigg(r^{*2}\frac{\partial T^{*}}{\partial r^{*}}\Bigg) - \frac{\partial T^{*}}{\partial t^{*}} - \frac{R_0^2\omega\rho_b C_{pb}}{k}T^{*} = -\frac{\dot{Q}R_0^2}{kT_{\infty}} - \frac{\dot{q_m}R_0^2}{kT_{\infty}} \end{align} \)

By defining the constants \(Q_0 = \frac{kT_{\infty}}{R_0^2}\) and \(C_0 = \frac{k}{R_0^2}\) that dimensionalize \(\dot{Q}\), \(\dot{q_m}\), and \(\omega\rho_b C_{pb}\) as \(Q^* = \frac{\dot{Q} + \dot{q_m}}{Q_0}\) and \(c^{*} = \frac{\omega\rho_b C_{pb}}{C_0}\), the bioheat eqation becomes

\(\begin{align} \large \frac{1}{r^{*2}}\frac{\partial}{\partial r^{*}}\Bigg(r^{*2}\frac{\partial T^{*}}{\partial r^{*}}\Bigg) - \frac{\partial T^{*}}{\partial t^{*}} - c^{*}T^{*} = -Q^{*}(r^{*}) \end{align} \)

We will use the following initial condition: \(T^{*}(r^*, t=0) = 0\) and the following boundary conditions: the temperature must remain finit at the center \(T^{*}(0, t^{*}) = finite\) and the temperature at the external spherical surface will be maintained at \(T_{\infty}\) so that the temperature \(T^{*}(r^{*} = 1,t) = 0\).

Analytical Solution

To obtain an analytical solution, let split the dimensionless heat generation term \(Q^{*}(r)\) into a constant heat generation due to metabolism \(q_m^{*}\) and a constant heat source of radius \(r_1\) at the center \(q_p^{c*}(r)\):

\(\begin{align} Q_{(r)}^{*} = q_m^{*} + q_p^{c*}(r) \end{align} \)

Substituting this into our non-dimensionalized bioheat equation, we get

\(\begin{align} \large \frac{1}{r^{*2}}\frac{\partial}{\partial r^{*}}\Bigg(r^{*2}\frac{\partial T^{*}}{\partial r^{*}}\Bigg) - \frac{\partial T^{*}}{\partial t^{*}} - c^{*}T^{*} = -q_m^{*} - q_p^{c*}(r) \end{align} \)

Now, let \(T^{*} = T_1(r, t) + T_2(r)\) (and we’ll drop the * for ease of notation) to get:

\(\begin{align} \large \frac{1}{r^{2}}\frac{\partial}{\partial r}\Bigg(r^{2}\frac{\partial T_{1}}{\partial r}\Bigg) - \frac{\partial T_1}{\partial t} - c T_1 = -q_p^{c}(r) \end{align} \)

\(\begin{align} \large \frac{1}{r^{2}}\frac{\partial}{\partial r}\Bigg(r^{2}\frac{\partial T_{2}}{\partial r}\Bigg) - c T_2 = -q_m \end{align} \)

Here are our boundary conditions for the above equations:

\(\begin{align} \large T_1(r=1, t) = 0, T_1(r=0,t) = finite \end{align} \)

\(\begin{align} \large T_2(r=1) = 0, T_2(r=0,t) = finite \end{align} \)

And the following is our initial condition:

\(\begin{align} \large T_1(r, 0) = -T_2(r) since T^{*}(r, 0) = T_1(r, 0) + T_2(r) = 0 \end{align} \)

Now, we’ll let \(T_1(r, t) = e^{-ct}U(r, t)\) and substitute into equation 8 above to get:

\(\begin{align} \large \frac{1}{r^{2}}\frac{\partial}{\partial r}\Bigg(r^{2}\frac{\partial U}{\partial r}\Bigg) - \frac{\partial U}{\partial t} = -q_p^{c}(r)e^{ct} = g(r, t) \end{align} \)

Now our source term is a function of \(r\) and \(t\), giving us \(g(r, t) = -q_p^c(r)e^{ct}\). The boundary conditions in terms of U becomes:

\(\begin{align} \large U(r = 1, t) = 0, U(r = 0, t) = finite \end{align} \)

And the initial conditions become:

\(\begin{align} \large U(r, t = 0) = -T_2(r) \end{align} \)

With a change of variable to \(T_2(r) = \frac{H(r)}{\sqrt{r}}\),

our bioheat equation 10 becomes a modified Bessle Equation in terms of H and the solution is obtained in terms of the modified Bessel Functions

\(\begin{align} \large T_2(r) = \frac{q_m}{c}\Bigg(1- \frac{I_{1/2}(r\sqrt{c})}{\sqrt{r}I_{1/2}(\sqrt{c})}\Bigg) \end{align} \)

The equation for U is obtained using Greens functions with an initial condition given by

\(\begin{align} \large U(r, 0) = -T_2(r) = -\frac{q_m}{c}\Bigg(1- \frac{I_{1/2}(r\sqrt{c})}{\sqrt{r}I_{1/2}(\sqrt{c})}\Bigg) \end{align} \)

A spherical heat source (eg, a magnetic nanoparticle solution) of radius \(r_1\) is placed at the cneter of the spherical tissue model and is represented as an impulse function with:

\(\begin{align} \large g(r, t) = g_p^c e^{ct} \frac{1}{4\pi r^2} \delta(r - r_1) \end{align} \)

Our solution for U will be found through

\(\begin{align} \large U(r, t) = \dfrac {2} {r}\sum _{m=1}^{\infty }e^{-\beta _{m}^{2}t}\sin \left( \beta _{m}r\right) \int _{0}^{1}r'\sin \left( \beta _{m}r'\right) \left( -T_{2}\left( r'\right) \right) dr' + \dfrac {1} {2\pi r_{1}r}\sum _{m=1}^{\infty }e^{-\beta _{m}^{2}t}\sin \left( \beta_m r\right)\sin \left( \beta _{m}r\right) \int _{\tau =0}^{t}e^{\beta _{m}^{2}\tau }g_{p}^{c}e^{c\tau }d\tau \end{align} \)

The integral \(\int _{\tau =0}^{t}e^{\beta _{m}^{2}\tau }g_{p}^{c}e^{c\tau }d\tau\) can be solved analytically for a constant intensity heat generation as

\(\begin{align} \large \int _{\tau =0}^{t}e^{\beta _{m}^{2}\tau }g_{p}^{c}e^{c\tau }d\tau = \dfrac {e^{\left( \beta _{m}^{2}+c\right) t}-1} {\beta _{m}^{2}+c}g_{p}^{c} \end{align} \)

With \(T_1 = Ue^{-ct}\)