# A Semi-Analytical Heterogeneous Model for Thermal Analysis of Cancerous Breasts

A. Ramírez-Torres, Reinaldo Rodríguez-Ramos, Aura Conci, Federico J.

Sabina, C. García-Reimbert, Luigi Preziosi, J Merodio, Frédéric Lebon

► To cite this version:

A. Ramírez-Torres, Reinaldo Rodríguez-Ramos, Aura Conci, Federico J. Sabina, C. García-Reimbert, et al.. A Semi-Analytical Heterogeneous Model for Thermal Analysis of Cancerous Breasts. Eddie Y. K. Ng; Mahnaz Etehadtavakol. Application of Infrared to Biomedical Sciences, Springer, pp.175-190, 2017, Series in BioEngineering, 10.1007/978-981-10-3147-2_11. hal-01697070

HAL Id: hal-01697070

https://hal.science/hal-01697070

Submitted on 3 Sep 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
# A Semi-Analytical Heterogeneous Model for Thermal Analysis of Cancerous Breasts

A. Ramírez-Torres, R. Rodríguez-Ramos, A. Conci, F.J. Sabina, C. García-Reimbert, L. Preziosi, J. Merodio and F. Lebon

**Abstract** In the present work coupled stationary bioheat transfer equations are considered. The cancerous breast is characterized by two areas of dissimilar thermal properties: the glandular and tumor tissues. The tumorous region is modeled as a two-phase composite where parallel periodic isotropic circular fibers are embedded in the glandular isotropic matrix. The periodic cell is assumed square. The local problem on the periodic cell and the homogenized equation are stated and solved. The temperature distribution of the cancerous breast is found through a numerical computation. A mathematical and computational model is integrated by FreeFem++.

**Keywords** Temperature · Cancerous breast · Asymptotic homogenization · Finite element method

A. Ramírez-Torres (✉) · F.J. Sabina · C. García-Reimbert

Instituto de Investigaciones en Matemáticas Aplicadas y en Sistemas, Universidad Nacional Autónoma de México, CP 20-126, 01000 CDMX, Mexico
e-mail: ariel.ramirez@iimas.unam.mx

F.J. Sabina

e-mail: fjs@mym.iimas.unam.mx

C. García-Reimbert

e-mail: cgr@mym.iimas.unam.mx

R. Rodríguez-Ramos

Facultad de Matemática y Computación, Universidad de La Habana,
CP 10400 La Habana, Cuba
e-mail: reinaldo@matcom.uh.cu

A. Conci

Instituto de Computação, Universidade Federal Fluminense,
CEP:24210-346 Rio de Janeiro, Brazil
e-mail: aconci@ic.uff.br

L. Preziosi

Dipartimento di Matematica, Politecnico di Torino, CP 10129 Turin, Italy
e-mail: luigi.preziosi@polito.it
# 1 Introduction

Actually, clinical examination, ultrasound, mammography, thermography, among others, are employed to identify and treat breast cancer [1, 2]. In particular, mammography is considered the standard procedure for detecting breast cancer. Yet, it presents difficulties for finding tumors in dense breasts. Thermography technique has arisen as a prospective method with the aim of increase the efficacy of the early discovery of breast cancer [3, 4]. Then, mathematical and numerical models have been proposed for studying thermal distribution on healthy and cancerous breasts, with the aim of using thermography as a complementary tool. For instance, [5] modeled a three-dimensional tumorous breast and sensitivity parameters are analyzed. Moreover, [3] were able to set a method to approximate thermal properties, where the physical process was ruled by a bioheat transfer equation. A three-dimensional breast, taking into account thermal and elastic properties, was modeled and the influence of both properties on the surface temperature was considered by [6]. In the aforementioned works, the numerical simulation was performed via FEM.

In the present study, a semi-analytical method is proposed for studying the breast thermal properties for different parameter data. Then, mathematical and computational modeling are integrated for solving two coupled stationary bioheat transfer equations. To separate micro and macro variables of the heterogeneous problem, the two-scale asymptotic expansion is used [7, 8]. In fact, multiscales methods have been successfully applied to various physical systems. For example, a formal two-scale asymptotic expansion for studying the macroscopic behavior of a porous and linear elastic solid was used in [9]. On the other hand, the homogeneous problem associated with the healthy breast tissues (without tumor) and the homogenized problem resulting by the application of the two-scale homogenization method to the heterogeneous tumor tissue, are solved using *FreeFem++*. Finally, numerical results are shown and discussed.

J. Merodio

Departamento de Mecánica de los Medios Continuos y T. Estructuras, E.T.S. de Caminos, Canales y Puertos, Universidad Politécnica de Madrid, CP 28040 Madrid, Spain
e-mail: jose.merodio@upm.es

F. Lebon

Laboratoire de Mécanique et d'Acoustique, Université Aix-Marseille, CNRS, Centrale Marseille, 4 Impasse Nikola Tesla, CS 40006, 13453 Marseille Cedex 13, France
e-mail: lebon@lma.cnrs-mrs.fr
## 2 Mathematical Model

The aim of the present work is to find the stationary temperature fields $u$ and $u^e$ that are described by the following bioheat transfer equations [10].

$$ egin{cases} -\frac{\partial}{\partial x_i} \left( K_{ij}^g \frac{\partial u}{\partial x_j} \right) + \rho_b c_b \omega_b^g u = q_m^g + \rho_b c_b \omega_b^g u_a & \text{in } \Omega_1, \\ -K_{ij}^g \frac{\partial u}{\partial x_j} n_i = h(u - u_e) & \text{on } \partial \Omega_1^n, \\ u = u_c & \text{on } \partial \Omega_1^d. \end{cases} \quad (1) $$

$$ egin{cases} -\frac{\partial}{\partial x_i} \left( K_{ij}^e \frac{\partial u^e}{\partial x_j} \right) + \rho_b c_b \omega_b^e u^e = q_m^e + \rho_b c_b \omega_b^e u_a & \text{in } \Omega_2, \\ K_{ij}^e \frac{\partial u^e}{\partial x_j} n_i = -K_{ij}^g \frac{\partial u_g}{\partial x_j} n_i & \text{on } \partial \Omega_2, \\ u^e = u_g & \text{on } \partial \Omega_2. \end{cases} \quad (2) $$

where $K_{ij}^g = k_g \delta_{ij}$ denotes glandular tissue thermal conductivity, $\rho_b$ is blood mass density, $c_b$ blood specific heat capacity, $u_a$ is the arterial blood temperature, $u_c$ the temperature at the boundary between breast and chest, $u_e$ is the surrounding temperature and $h$ represents the combined effective heat transfer coefficient due to convection, radiation and evaporation of 13.5 W/m² K [11]. Besides, the rapidly oscillating coefficients $K_{ij}^e, \omega_b^e$ and $q_m^e$ are defined as follows

$$ K_{ij}^{e}(\mathbf{x}) = \begin{cases} k_g \delta_{ij}, & \mathbf{x} \in \Omega_g^e, \\ k_t \delta_{ij}, & \mathbf{x} \in \Omega_t^e, \end{cases} \quad \omega_b^e(\mathbf{x}) = \begin{cases} \omega_b^g, & \mathbf{x} \in \Omega_g^e, \\ \omega_b^t, & \mathbf{x} \in \Omega_t^e \end{cases}, \quad \text{and} \quad q_m^e(\mathbf{x}) = \begin{cases} q_m^g, & \mathbf{x} \in \Omega_g^e, \\ q_m^t, & \mathbf{x} \in \Omega_t^e \end{cases}. $$

Note that in the case of the healthy breast model only ($P_1$) has to be solved.

For the sake of simplicity, we will work in a two-dimensional section where the breast geometry is represented by a hemispherical shape with a diameter $L$ as done in [11]. The healthy breast will be represented by a homogeneous tissue (glandular tissue) and associated with the open, bounded, and connected domain $\Omega_1$ with Lipschitz boundary $\partial \Omega_1 = \partial \Omega_1^n \cup \partial \Omega_1^d$, where $\partial \Omega_1^n \cap \partial \Omega_1^d = \emptyset$. On the other hand, the cancerous tissue will be characterized by two regions of dissimilar thermal properties: the tumoral area ($\Omega_t^e$—fibers) and the glandular area ($\Omega_g^e$—matrix). In this sense, the cancerous region will consist of a periodic microstructure associated with the open, bounded, and connected domain $\Omega_2 = \Omega_g^e \cup \Omega_t^e \cup \partial \Omega_g^e$ with Lipschitz boundary $\partial \Omega_2 = \partial \Omega_g^e$ and with $\Omega_g^e \cap \Omega_t^e = \emptyset$. Then, the cancerous breast is represented by $\Omega = \Omega_2 \cup \Omega_1$ (Fig. 1). Let $\varepsilon > 0$ be the size of the microstructure and $\mathbf{y} = \mathbf{x}/\varepsilon$ the fast scale coordinate. The reference periodic cell will be denoted by $Y$, which contains one inclusion occupying the domain $Y_t$ with Lipschitz boundary $\partial Y_t$ such that $Y = Y_g \cup Y_t \cup \partial Y_t$, with $\bar{Y}_t \subset Y$ and $Y_g \cap Y_t = \emptyset$. It is also assumed that $\Omega_g^e$ is connected and that the inclusions do not intersect the boundary $\partial \Omega_g^e$. In previous works, soft tissues assume to present this type of arrangement. In fact, Penta et al. [9]
<figure><img src="image_2.png" /><figcaption>**Fig. 1** Decomposition of the macroscopic domain (**left**) and the corresponding unit periodic cell (**right**)</figcaption></figure>

used the same periodic geometry to depict a porous tissue microstructure. Boundary conditions for (1) are heat transfer by convection between the surface of the tissue and the external environment on $\partial\Omega_1^n$ and a prescribed temperature on $\partial\Omega_1^d$. In the case of (2) we assume heat and temperature continuity on $\partial\Omega_2$. Moreover, continuity conditions for temperature and heat flow are imposed on $\Gamma^\epsilon$ (boundary between the glandular tissue $\Omega_g^\epsilon$ and the tumor inclusions $\Omega_t^\epsilon$), i.e.,

$$ \|u^\epsilon\| = 0 \quad \text{on } \Gamma^\epsilon \quad \text{and} \quad \|\mathbf{K}^\epsilon \nabla_x u^\epsilon \cdot \mathbf{n}\| = 0 \quad \text{on } \Gamma^\epsilon. \quad (3) $$

### 3 Two-Scale Homogenization

Here, the two-scale homogenization technique is applied to find the homogenized equation and corresponding effective coefficients. An overview of how this method is applied and its main assumptions can be found in [12]. Specifically, after finding the solution $u$ of problem (1), an asymptotic expansion of $u^\epsilon$ [problem (2)] is sought as a function of $\epsilon$ for $\epsilon \to 0$, namely

$$ u^{\epsilon}(\mathbf{x}) = u^{0}(\mathbf{x}) + \epsilon\chi_{p}(\mathbf{y}) \frac{\partial u^{0}(\mathbf{x})}{\partial x_{p}} + \epsilon^{2}u^{2}(\mathbf{x}, \mathbf{y}) + \dots, \quad (4) $$

where the functions $\chi_p(\mathbf{y})$, $u^2(\mathbf{x}, \mathbf{y})$, are $Y$-periodic in $\mathbf{y}$. In particular, the vector function $\chi(\mathbf{y})$ satisfies the unit cell problem

$$ \left\{ \begin{array}{ll} -\frac{\partial}{\partial y_i} \left( K_{ij}(\mathbf{y}) \frac{\partial \chi_p(\mathbf{y})}{\partial y_j} + K_{ip}(\mathbf{y}) \right) = 0 & \text{in } Y \setminus \Gamma, \\ \|\chi\| = 0 & \text{on } \Gamma, \\ \left[ \left( K_{ij}(\mathbf{y}) \frac{\partial \chi_p(\mathbf{y})}{\partial y_j} + K_{ip}(\mathbf{y}) \right) n_i \right] = 0 & \text{on } \Gamma \end{array} \right. \quad (5) $$

and $u^0(\mathbf{x})$ the homogenized problem solution
$$
\left\{
\begin{array}{ll}
-\frac{\partial}{\partial x_i} \left( \hat{K}_{ij} \frac{\partial u^0(\mathbf{x})}{\partial x_j} \right) + \left( \rho_b c_b \omega_b^g \frac{|Y_g|}{|Y|} + \rho_b c_b \omega_b^t \frac{|Y_t|}{|Y|} \right) u^0(\mathbf{x}) = (q_m^g + \rho_b c_b \omega_b^g u_a) \frac{|Y_g|}{|Y|} + & \text{in } \Omega_2, \\
+ (q_m^t + \rho_b c_b \omega_b^t u_a) \frac{|Y_t|}{|Y|} & \text{on } \partial\Omega_2, \\
u^0(\mathbf{x}) = u_g &
\end{array}
\right.
\tag{6}
$$

where $|\cdot|$ represents volume fraction. The effective constant coefficients $\hat{K}_{ip}$ are
given by

$$
\hat{K}_{ip} = \left\langle K_{ij} \frac{\partial \chi_p}{\partial y_j} + K_{ip} \right\rangle \quad (7)
$$

where $p = 1, 2$ and $\langle \cdot \rangle$ denotes volume average.

**3.1 Analytical Solution of the Unit Cell Problem**

In particular, the theory of analytical functions by Muskhelishvili [13] is applied to
solve the cell problem (5). In this sense, the solutions of the local problems are
written as

$$
\chi_1^{(g)} = \operatorname{Re} \left\{ a_0^1 z + \sum_{k=1}^{\infty} a_k^1 \frac{\zeta^{(k-1)}(z)}{(k-1)!} \right\} \quad \text{and} \quad \chi_1^{(t)} = \operatorname{Re} \left\{ \sum_{l=1}^{\infty} c_l^1 z^l \right\}, \quad (8)
$$

$$
\chi_2^{(g)} = \operatorname{Im} \left\{ a_0^2 z + \sum_{k=1}^{\infty} a_k^2 \frac{\zeta^{(k-1)}(z)}{(k-1)!} \right\} \quad \text{and} \quad \chi_2^{(t)} = \operatorname{Im} \left\{ \sum_{l=1}^{\infty} c_l^2 z^l \right\}, \quad (9)
$$

where $(\gamma)$ with $\gamma = g, t$ denotes the constituent, the superscript $o$ specifies that the
sum is carried out over odd indices, the coefficients $a_0^p, a_k^p$ and $c_l^p (p = 1, 2)$ are real
and $\zeta$ is the zeta quasi periodic Weierstrass function. Now, using Laurent's
expansion of $\zeta$ and the quasi-periodicity property of $\zeta$ and its derivatives

$$
\chi_1^{(g)} = \operatorname{Re} \left\{ \sum_{l=1}^{\infty} (a_l^1 z^{-l} - A_l^1 z^l) \right\} \quad \text{and} \quad \chi_2^{(g)} = \operatorname{Im} \left\{ \sum_{l=1}^{\infty} (a_l^2 z^{-l} - A_l^2 z^l) \right\}, \quad (10)
$$

where for $p = 1, 2$

$$
A_l^p = \sum_{k=1}^{\infty} k a_k^p \eta_{kl} \quad \text{with} \quad \eta_{kl} = \begin{cases} (-1)^{p+1} \pi, & k+l=2 \\ \frac{(k+l-1)!}{k!l!} S_{k+l}, & k+l>2 \end{cases}
$$
and $S_k$ are called the reticulate sums and are defined as $S_k = \sum_{w \in L^*} \frac{1}{w^k} (k \ge 3, k$ odd) with $L^*$ representing the lattice excluding the number $w = 0$ and $w = mw_1 + nw_2$ where $m, n \in \mathbb{Z}$ and $w_1, w_2$ are the periods. In particular, in the present work $w_1 = 1$ and $w_2 = i$, due to we are in presence of square unit cells.

Substitution of (8)–(10) in boundary conditions of problem (5) and taking into account that on $\Gamma, z = Re^{i\theta}$ where $R$ is the circumference radius give that coefficients $a_k^p$ can be found through solution of the following infinite linear system (for finding the effective properties it is truncated into an appropriate order $k = N$)

$$ (\xi^{-1}\mathbf{I} + (-1)^{p+1}\mathbf{W}^p)\hat{\mathbf{A}}^p = \mathbf{V}^p, \qquad (11) $$

where $\hat{\mathbf{A}}^p = (\hat{a}_1^p, \hat{a}_2^p, \dots)^T$, $\hat{a}_k^p = a_k^p \sqrt{k}/R^k$, $\mathbf{V}^p = ((-1)^{p+1}R, 0, \dots)^T$,

$$ \xi^{-1} = \frac{k_g + k_t}{k_g - k_t} \quad \text{and} \quad \mathbf{W}^p = \begin{cases} (-1)^{p+1} \pi R^2, & k+l=2 \\ \sum_{k=1}^{\infty} \sqrt{k} \sqrt{l} \eta_{kl} R^{k+l}, & k+l>2. \end{cases} $$

Now, from Eq. (7)

$$ \hat{K}_{ip} = \left\langle K_{ij} \frac{\partial \chi_p}{\partial y_j} + K_{ip} \right\rangle. $$

Using the form of $K_{ij}$, Green's theorem, the double periodicity of $\chi_p$ and formulas (8)–(10), then

$$ \hat{K}_{pp} = \begin{cases} k_g(1 - 2\pi a_1^2), & \text{if } p = 1 \\ k_t(1 + 2\pi a_1^2), & \text{if } p = 2. \end{cases} \qquad (12) $$

In fact, if $k_g = k_t$. Then, $\hat{K} = \hat{K}_{11} = \hat{K}_{22}$.

# 4 Numerical Solution and Analysis of Results

This section is devoted to find the temperature distribution of problems (1) and (6) where we define as $g = \rho_b c_b \omega_b^g$, $g^\epsilon(\mathbf{x}) = \rho_b c_b \omega_b^g(\mathbf{x})$, $f = q_m^g + \rho_b c_b \omega_b^g u_a$, $f^\epsilon(\mathbf{x}) = q_m^\epsilon(\mathbf{x}) + \rho_b c_b \omega_b^\epsilon(\mathbf{x}) u_a$. With this aim, we follow the following procedure.
## 4.1 Solution of (1)

For finding solution $u$ of (1), we use FreeFem++. First, the problem must be written in its weak formulation. In this sense, let $H_d^1(\Omega_1) = \{u \in H^1(\Omega_1) \text{ s.t. } \gamma(u) = 0 \text{ on } \partial\Omega_1^d\}$. Using the trace theorem for $u_c \in H^{1/2}(\Omega_1)$, there exists a continuous linear operator $R_0: H_\frac{1}{2}(\partial\Omega_1) \to H^1(\Omega_1)$ such that $\gamma(R_0u_c) = u_c$. Now, we define $\tilde{u} = u - R_0u_c \in H_d^1(\Omega_1)$. Then, on $\partial\Omega_1^d$

$$ \gamma(\tilde{u}) = \gamma(u) - \gamma(R_0u_c) = u_c - u_c = 0. $$

In this way, the equivalent variational formulation of problem (1) is

$$ \left\{ \begin{array}{ll} \text{Find } \tilde{u} \in H_d^1(\Omega_1) \text{ such that} & \\ a(\tilde{u}, v) = L(v), \quad \forall v \in H_d^1(\Omega_1), & \end{array} \right. \qquad (13) $$

where

$$ a(\tilde{u}, v) = \int_{\Omega_1} \mathbf{K}^g \nabla_x \tilde{u} \cdot \nabla_x v \, dx + \int_{\Omega_1} g \tilde{u} v \, dx + \int_{\partial \Omega_1^\eta} h \tilde{u} v \, dS, $$

$$ L(v) = \int_{\Omega_1} fv \, dx - \int_{\Omega_1} g(R_0 u_c) v \, dx - \int_{\Omega_1} \mathbf{K}^g \nabla_x (R_0 u_c) \cdot \nabla_x v \, dx \\ + \int_{\partial \Omega_1^\eta} h(u_c - R_0 u_c) v \, dS. $$

In particular, the weak solution existence and uniqueness of problem (13) can be proved by standard methods using the Lax–Milgram theorem. In this sense, the following must be proved:

(i) The bilinear form $a(\tilde{u}, v)$ is continuous

In this sense, observe that $\mathbf{K}^g \in M(\alpha, \beta, \Omega_1)$ and by Cauchy–Schwartz

$$ |a(\tilde{u}, v)| \le \beta \|\nabla_x \tilde{u}\|_{L^2(\Omega_1)} \|\nabla_x v\|_{L^2(\Omega_1)} + g \|\tilde{u}\|_{L^2(\Omega_1)} \|v\|_{L^2(\Omega_1)} \\ + h \|\tilde{u}\|_{L^2(\partial \Omega_1^\eta)} \|v\|_{L^2(\partial \Omega_1^\eta)}. $$

Furthermore, by the Poincaré–Friedrichs II theorem

$$ \|\tilde{u}\|_{L^2(\Omega_1)} \le \tilde{C}_1 \|\nabla_x \tilde{u}\|_{L^2(\Omega_1)} = \tilde{C}_1 \|\tilde{u}\|_{H_d^1(\Omega_1)}. $$

On the other hand, by the trace theorem
$$
\begin{aligned}
\|\tilde{u}\|_{L^2(\partial\Omega_1^n)} &= \|\gamma(\tilde{u})\|_{L^2(\partial\Omega_1^n)} = \left[ \int_{\partial\Omega_1^n} |\gamma(\tilde{u})|^2 + \int_{\partial\Omega_1^d} |\gamma(\tilde{u})|^2 \right]^{\frac{1}{2}} \\
&= \|\tilde{u}\|_{L^2(\partial\Omega_1)} \le C_2 \|\tilde{u}\|_{H^1(\Omega_1)}.
\end{aligned}
$$

Moreover,

$$
\begin{aligned}
\|\tilde{u}\|_{H^1(\Omega_1)} &= \left[ \|\tilde{u}\|_{L^2(\Omega_1)}^2 + \|\nabla_x \tilde{u}\|_{L^2(\Omega_1)}^2 \right]^{\frac{1}{2}} \le \left[ \tilde{C}_1 \|\nabla_x \tilde{u}\|_{L^2(\Omega_1)}^2 + \|\nabla_x \tilde{u}\|_{L^2(\Omega_1)}^2 \right]^{\frac{1}{2}} \\
&= (1 + \tilde{C}_1)^{\frac{1}{2}} \|\nabla_x \tilde{u}\|_{L^2(\Omega_1)} = \tilde{C}_3 \|\tilde{u}\|_{H_d^1(\Omega_1)}, \quad \tilde{u} \in H_d^1(\Omega_1).
\end{aligned}
$$

Therefore,

$$
\| \tilde{u} \|_{L^2(\partial \Omega_1^n)} \le C_2 \| \tilde{u} \|_{H^1(\Omega_1)} \le \tilde{C}_4 \| \tilde{u} \|_{H_d^1(\Omega_1)},
$$

where $\tilde{C}_4 = C_2 \tilde{C}_3 = C_2(1 + \tilde{C}_1)^{\frac{1}{2}}$. Finally, for $\tilde{u} \in H_d^1(\Omega_1)$

$$
\begin{aligned}
|a(\tilde{u}, v)| &\le \beta \|\tilde{u}\|_{H_d^1(\Omega_1)} \|v\|_{H_d^1(\Omega_1)} + g \tilde{C}_1 C_1 \|\tilde{u}\|_{H_d^1(\Omega_1)} \|v\|_{H_d^1(\Omega_1)} + h \tilde{C}_4 C_4 \|\tilde{u}\|_{H_d^1(\Omega_1)} \|v\|_{H_d^1(\Omega_1)} \\
&\le C_5 \|\tilde{u}\|_{H_d^1(\Omega_1)} \|v\|_{H_d^1(\Omega_1)},
\end{aligned}
$$

with $C_5 = \beta + g\tilde{C}_1C_1 + h\tilde{C}_4C_4$.

(ii) The bilinear form $a(\tilde{u}, v)$ is $H_d^1$-elliptic

Let $u \in H_d^1(\Omega_1)$,

$$
\begin{aligned}
a(u, u) &= \int_{\Omega_1} \mathbf{K}^g (\nabla_x u)^2 dx + \int_{\Omega_1} gu^2 dx + \int_{\partial \Omega_1^n} hu^2 dS \\
&\geq C_6 \left( \int_{\Omega_1} (\nabla_x u)^2 dx + \int_{\Omega_1} u^2 dx + \int_{\partial \Omega_1^n} u^2 dS \right) \\
&= C_6 \left( \|\nabla_x u\|_{L^2(\Omega_1)}^2 + \|u\|_{L^2(\Omega_1)}^2 + \|u\|_{L^2(\partial \Omega_1^n)}^2 \right) \\
&\geq C_6 \|\nabla_x u\|_{L^2(\Omega_1)}^2 = C_6 \|u\|_{H_d^1(\Omega_1)}^2, \quad \text{with } C_6 = \min(\alpha, g, h).
\end{aligned}
$$

(iii) The linear form $L(v)$ is continuous in $H_d^1(\Omega_1)$

Let $v \in H_d^1(\Omega_1)$,
$$
\begin{aligned}
|L(v)| &= \left| \int_{\Omega_1} f v \, dx - \int_{\Omega_1} g(R_0 u_c) v \, dx - \int_{\Omega_1} \mathbf{K}^g \nabla_x (R_0 u_c) \cdot \nabla_x v \, dx + \int_{\partial \Omega_1^n} h(u_e - R_0 u_c) v \, dS \right| \\
&\le \int_{\Omega_1} |f v| \, dx + \int_{\Omega_1} |g(R_0 u_c) v| \, dx + \int_{\Omega_1} |\mathbf{K}^g \nabla_x (R_0 u_c) \cdot \nabla_x v| \, dx + \int_{\partial \Omega_1^n} |h(u_e - R_0 u_c) v| \, dS.
\end{aligned}
$$

But, by Cauchy–Schwarz inequality, the Poincaré–Friedrichs II theorem and the fact that $R_0 u_c \in H^1(\Omega_1)$ and $\nabla_x(R_0 u_c) \in (L^2(\Omega_1))^n$

$$|L(v)| \le C_{10} \|v\|_{H_d^1(\Omega_1)}$$

where $C_{10} = C_1 \|f\|_{L^2(\Omega_1)} + gC_1 \|R_0 u_c\|_{L^2(\Omega_1)} + \beta \|\nabla_x(R_0 u_c)\|_{L^2(\Omega_1)} + (hC_7|u_e| + hC_8\|u_c\|_{H^{1/2}(\partial\Omega_1^n)})C_3$.

Thus, (i)–(iii) proves the existence and uniqueness of solution $\tilde{u}^0$ by using the Lax–Milgram theorem.

Now, it must be shown that the map $\langle f^e \rangle \in L^2(\Omega_1) \to u \in H_d^1(\Omega_1)$ is continuous in order to prove the regularity of the weak solution. In fact, from the $H_d^1$-ellipticity of the bilinear form

$$|a(u, u)| \ge C_6 \|u\|_{H_d^1(\Omega_1)}^2$$

and the continuity of the linear operator in $H_d^1(\Omega_1)$

$$|L(u)| \le C_{10} \|u\|_{H_d^1(\Omega_1)}^2.$$

Then,

$$C_6 \|u\|_{H_d^1(\Omega_1)}^2 \le |a(u, u)| = |L(u)| \le C_{10} \|u\|_{H_d^1(\Omega_1)}^2,$$

i.e.

$$\|u\|_{H_d^1(\Omega_1)} \le \frac{C_{10}}{C_6} = C_{11} \|f\|_{L^2(\Omega_1)}.$$

with

$$C_{10} = \frac{1}{C_6} \left( C_1 + \frac{g C_1 \|R_0 u_c\|_{L^2(\Omega_1)}}{\|f\|_{L^2(\Omega_1)}} + \frac{\beta \|\nabla_x(R_0 u_c)\|_{L^2(\Omega_1)}}{\|f\|_{L^2(\Omega_1)}} + \frac{(hC_7|u_e| + hC_8\|u_c\|_{H^{1/2}(\partial\Omega_1^n)})C_3}{\|f\|_{L^2(\Omega_1)}} \right).$$

Now, the contribution of $R_0 u_c$ may be difficult in some cases. However, FreeFem++ replaces the Dirichlet condition by a Robin condition of the form
$\nabla_x u \cdot n + u/\epsilon = u_c/\epsilon$ on $\partial\Omega_1^d$ and solves the problem with a very small value of $\epsilon$. In particular, we approximate the involved functions by piecewise linear continuous finite elements.

## 4.2 Solution of (6)

The last step in the homogenization procedure is to solve the homogenized problem (6). Here we prove that $u^0$ is its solution and that the problem is well posed. In this sense, let $H_0^1(\Omega_2) = \{u \in H^1(\Omega_2) \text{ s.t. } \gamma(u) = 0 \text{ on } \partial\Omega_2\}$. Using the trace theorem for $u_g \in H^{1/2}(\Omega_2)$ there exists a continuous linear operator $R_0: H^1(\partial\Omega_2) \to H^1(\Omega_2)$ such that $\gamma(R_0u_g) = u_g$. Now, we define $\tilde{u}^0 = u^0 - R_0u_g \in H_0^1(\Omega_2)$. Then, on $\partial\Omega_2$

$$ \gamma(\tilde{u}^0) = \gamma(u^0) - \gamma(R_0 u_g) = u_g - u_g = 0. $$

In this way, the equivalent variational formulation of problem (6) is

$$ \left\{ \begin{array}{l} \text{Find } \tilde{u}^0 \in H_0^1(\Omega_2) \text{ such that} \\ a(\tilde{u}^0, v) = L(v), \quad \forall v \in H_0^1(\Omega_2) \end{array} \right. , \qquad (14) $$

where

$$ \begin{aligned} a(\tilde{u}^0, v) &= \int_{\Omega_2} \hat{\mathbf{K}} \nabla_x \tilde{u}^0 \cdot \nabla_x v \,dx + \int_{\Omega_2} \langle g^\epsilon \rangle \tilde{u}^0 v \,dx, \\ L(v) &= \int_{\Omega_2} \langle f^\epsilon \rangle v \,dx - \int_{\Omega_2} \hat{\mathbf{K}} \nabla_x (R_0 u_g) \cdot \nabla_x v \,dx - \int_{\Omega_2} \langle g^\epsilon \rangle (R_0 u_g) v \,dx. \end{aligned} $$

In particular, $\hat{\mathbf{K}} \in M(\alpha, \beta, \Omega_2)$ see Cionarescu and Donato [14], and $\langle f^\epsilon \rangle \in L^2(\Omega_2)$ and $\langle g^\epsilon \rangle > 0$. The existence and uniqueness of solution $\tilde{u}^0$ can be proved through the Lax-Milgram theorem. Then, we must show that:

(i) The bilinear form $a(\tilde{u}^0, v)$ is continuous

In this sense, observe that using the fact that $\hat{\mathbf{K}} \in M(\alpha, \beta, \Omega_2)$ and by Cauchy-Schwartz inequality

$$ |a(\tilde{u}^0, v)| \le \beta \|\nabla_x \tilde{u}^0\|_{L^2(\Omega_2)} \|\nabla_x v\|_{L^2(\Omega_2)} + \langle g^\epsilon \rangle \|\tilde{u}^0\|_{L^2(\Omega_2)} \|v\|_{L^2(\Omega_2)}. $$

Now, from remark 3.37 p. 32 by Cionarescu and Donato [14], for $\tilde{u}^0, v \in H_0^1(\Omega_2)$,
$$ \|\tilde{u}^0\|_{H_0^1(\Omega_2)} = \|\nabla_x \tilde{u}^0\|_{L^2(\Omega_2)}, $$

$$ \|v\|_{H_0^1(\Omega_2)} = \|\nabla_x v\|_{L^2(\Omega_2)}. $$

Furthermore, by the Poincaré–Friedrichs I theorem,

$$ \|\tilde{u}^0\|_{L^2(\Omega_2)} \le \hat{C}_1 \|\nabla_x \tilde{u}^0\|_{L^2(\Omega_2)} = \hat{C}_1 \|\tilde{u}^0\|_{H_0^1(\Omega_2)}, $$

$$ \|v\|_{L^2(\Omega_2)} \le C_1 \|\nabla_x v\|_{L^2(\Omega_2)} = C_1 \|v\|_{H_0^1(\Omega_2)}. $$

Finally, for $\tilde{u}^0 \in H_0^1(\Omega_2)$

$$ |a(\tilde{u}^0, v)| \le \beta \|\tilde{u}^0\|_{H_0^1(\Omega_2)} \|v\|_{H_0^1(\Omega_2)} + \langle g^\epsilon \rangle \hat{C}_1 C_1 \|\tilde{u}^0\|_{H_0^1(\Omega_2)} \|v\|_{H_0^1(\Omega_2)} $$

$$ \le C_2 \|\tilde{u}^0\|_{H_0^1(\Omega_2)} \|v\|_{H_0^1(\Omega_2)}, $$

with $C_2 = \beta + \langle g^\epsilon \rangle \hat{C}_1 C_1$.

(ii) The bilinear form $a(\tilde{u}^0, v)$ is $H_0^1$-elliptic

Let $u \in H_0^1(\Omega_2)$,

$$ \begin{aligned} a(u, u) &= \int_{\Omega_2} \hat{\mathbf{K}}(\nabla_x u)^2 dx + \int_{\Omega_2} \langle g^\epsilon \rangle u^2 dx \\ &= C_3 \left( \|\nabla_x u\|_{L^2(\Omega_2)}^2 + \|u\|_{L^2(\Omega_2)}^2 \right), \quad \text{with } C_3 = \min(\alpha, \langle g^\epsilon \rangle) \\ &\ge C_3 \|\nabla_x u\|_{L^2(\Omega_2)}^2 = C_3 \|u\|_{H_0^1(\Omega_2)}^2. \end{aligned} $$

(iii) The linear form $L(v)$ is continuous in $H_0^1(\Omega_2)$

Let $v \in H_0^1(\Omega_2)$,

$$ |L(v)| \le \int_{\Omega_2} |\langle f^\epsilon \rangle v| dx + \int_{\Omega_2} |\hat{\mathbf{K}} \nabla_x (R_0 u_g) \cdot \nabla_x v| dx + \int_{\Omega_2} |\langle g^\epsilon \rangle (R_0 u_g) v| dx. $$

But, using the Cauchy–Schwarz inequality, the Poincaré–Friedrichs I theorem and the fact that $R_0 u_g \in H^1(\Omega_2)$ and $\nabla_x (R_0 u_g) \in (L^2(\Omega_2))^n$

$$ |L(v)| \le C_6 \|v\|_{H_0^1(\Omega_2)}, $$

where $C_6 = C_4 \|\langle f^\epsilon \rangle\|_{L^2(\Omega_2)} + \langle g^\epsilon \rangle C_5 \|R_0 u_g\|_{L^2(\Omega_2)} + \beta \|\nabla_x (R_0 u_g)\|_{L^2(\Omega_2)}$.

Thus, (i)–(iii) proves the existence and uniqueness of solution $\tilde{u}^0$ by using the Lax–Milgram theorem.
Now, it must be shown that the map $\langle f^\epsilon \rangle \in L^2(\Omega_2) \to u \in H_0^1(\Omega_2)$ is continuous in order to prove the regularity of the weak solution. In fact, from the $H_0^1$-ellipticity of the bilinear form

$$|a(u, u)| \geq C_3 \|u\|_{H_0^1(\Omega_2)}^2$$

and the continuity of the linear operator in $H_0^1(\Omega_2)$

$$|L(u)| \leq C_6 \|u\|_{H_0^1(\Omega_2)}^2.$$

Then,

$$C_3 \|u\|_{H_0^1(\Omega_2)}^2 \leq |a(u, u)| = |L(u)| \leq C_6 \|u\|_{H_0^1(\Omega_2)}^2,$$

i.e.,

$$\|u\|_{H_0^1(\Omega_2)} \leq \frac{C_6}{C_3} = C_7 \|\langle f^\epsilon \rangle\|_{L^2(\Omega_2)},$$

with

$$C_7 = \frac{1}{C_3} C_4 + \frac{\langle g^\epsilon \rangle C_5 \|R_0 u_g\|_{L^2(\Omega_2)}}{\|\langle f^\epsilon \rangle\|_{L^2(\Omega_2)}} + \frac{\beta \|\nabla_x (R_0 u_g)\|_{L^2(\Omega_2)}}{\|\langle f^\epsilon \rangle\|_{L^2(\Omega_2)}}.$$

Once the solution $u$ of (1) is found, we proceed to solve (6). In particular, the homogenized problem (6) is solved using the aforementioned *FreeFem++*. As above, we approximate the involved functions by piecewise linear continuous finite elements.

## 4.3 Analysis of Results

Numerical calculations are carried out for three breast models A, B, and C, whose tissue parameters are shown in Table 1. Temperatures are fixed as $u_a = u_c = 37^\circ C$ [5]. We fixed the surrounding temperature $u_e = 20$. The metabolic heat value for different tumor sizes follows the law given by Jiang et al. [6] as $q_m^t = C/(468.6 \ln(100D) + 50)$, where $C = 3.27 \times 10^6$ Wday/m$^3$ and $D$ is the tumor diameter.

Figure 2A–C show the temperature distribution of healthy breast tissues with $L = 0.14$ m, i.e., without a tumor. Now, in Fig. 3 it is shown how depth (in the present study the depth is stated as the distance between the tumor center and the point on the breast surface in the same axis) affects breast thermal distribution. In particular, in the zone “far” from the tumor area, no appreciable temperature
<table><caption>**Table 1** Tissue parameters</caption><thead><tr><th rowspan="2">Model</th><th colspan="2">k (W/m °C)</th><th colspan="2">ω<sub>b</sub> (1/s)</th><th rowspan="2">c<sub>b</sub> (J/Kg °C)</th><th rowspan="2">ρ<sub>b</sub> (kg/m<sup>3</sup>)</th><th rowspan="2">q<sub>g</sub> (Wm<sup>-3</sup>)</th><th rowspan="2">Reference</th></tr><tr><th>k<sub>t</sub></th><th>k<sub>g</sub></th><th>ω<sub>b</sub><sup>t</sup></th><th>ω<sub>b</sub><sup>g</sup></th></tr></thead><tbody><tr><td>A</td><td>0.48</td><td>0.48</td><td>0.0132</td><td>0.0006612</td><td>3300</td><td>1100</td><td>700</td><td>[5]</td></tr><tr><td>B</td><td>0.48</td><td>0.48</td><td>0.009</td><td>0.00018</td><td>4200</td><td>1060</td><td>450</td><td>[3]</td></tr><tr><td>C</td><td>0.511</td><td>0.48</td><td>0.0108</td><td>0.000539</td><td>4200</td><td>1060</td><td>700</td><td>[15]</td></tr></tbody></table>

<figure><img src="image_3.png" /><figcaption>**Fig. 2** Thermal distribution of a healthy breast tissue</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Fig. 3** Thermal distribution of a cancerous breast tissue with an embedded spherical tumor of radius $r = 0.01$ m located at depths $d = 0.03$ m (first line) and $d = 0.05$ m (second line), respectively</figcaption></figure>

changes at the surface are observed. Which is not the case when the tumor is located near to the boundary where the temperature difference at the surface is higher if compared with Fig. 2. Indeed, when the tumor is nearer to the boundary, the surface temperature increases. This behavior is in accordance with the observations made by [5] and [6]. Figure 4 first line of graphs, center line and bottom line, present the steady-state temperature for a cancerous breast tissue with $L = 0.13$ m, $L = 0.15$ m and $L = 0.17$ m, respectively. In particular, a sphere with radius $r = 0.01$ m was inserted in the breast model to imitate the in situ tumor at a depth of $d = 0.04$ m. Moreover, a relative large tumor volume fraction $|Y_t| = 0.7$ is considered so that healthy breast tissue volume fraction is $|Y_g| = 0.3$ in the tumorous region. As observed, if the breast dimensions are bigger, the maximum temperature is higher,
<figure><img src="image_5.png" /><figcaption>Fig. 4 Thermal distribution of a cancerous breast tissue with an embedded spherical tumor of radius $r = 0.01$ m for $L = 0.13$ m, $L = 0.15$ m and $L = 0.17$ m, respectively</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5 Thermal distribution of a cancerous breast tissue with an embedded spherical tumor of radius $r = 0.01$ m located $0.02$ m at the right and left of the central axis, respectively</figcaption></figure>

where the temperature decreases from the chest wall to the front breast. Moreover, surface temperature varies with breast dimension. On the other hand, the first and second line of Fig. 5 show the temperature for a cancerous breast tissue where the tumor is not located on the central axis $x_1 = 0.07$ m, i.e., its center is situated $0.02$ m at the right and left of the central axis, respectively. Even, when the tumor is found off central axis, it influences the temperature behavior on the nearest boundary, which is higher than that of the adjacent surface. Besides, for Figs. 4 and 5 a temperature variation is noticed in the tumor area and the region surrounding it.
# 5 Conclusions

Here, a semi-analytical method is used for studying breast thermography through coupled stationary bioheat transfer equations. One hand, the breast is assumed to be homogeneous and constituted by glandular tissue. On the other hand, the tumor area is represented by a periodic composite and comprised of glandular and cancerous tissue. In particular, the temperature distribution on both, breast and tumor tissue, was computed using a numerical algorithm implemented in FreeFem++. In summary, if the breast dimensions are bigger than the maximum temperature is higher and no appreciable changes in temperature difference were observed far from the breast boundary. The work results also indicate that the data parameter will influence the thermal distribution of the tumorous breast. The proposed method provides a helpful framework for studying the thermal profile of breast cancerous tissues. Moreover, it facilitates the understanding of the complex behavior of its surface temperature. Also, it improves the current premature discovery and analysis of breast tumors, integrating mathematical and computational tools. In fact, thermography together with mathematical and computational modeling bring an appropriate methodology in order to allow the assessment of rapidly growing neoplasm.

**Acknowledgements** AR gratefully acknowledges the Program of Postdoctoral Scholarships of DGAPA from UNAM, Mexico. The authors would like to thank Professor Sébastien Martin, Université Paris Descartes, Paris, France for helpful discussions during the course of this work. Also, we thank Departamento de Matemáticas y Mecánica IIMAS-UNAM for their support and Ramiro Chávez Tovar and Ana Pérez Arteaga for computational assistance. JM acknowledges support from the Ministerio de Ciencia in Spain under the project reference DPI2014-58885-R. Thanks to the Project (7515) Métodos Físico-Matemáticos para el estudio de nuevos materiales y la propagación de ondas. Aplicaciones. A.C. is grateful to CNPq (Projects INCT-MACC and 303240/2015-6)

# References

1. Muralidhar, G.S., Bovik, A.C., Sampat, M.P., Whitman, G.J., Haygood, T.M., Stephens, T.W., Markey, M.K.: Computer-aided diagnosis in breast magnetic resonance imaging. Mt Sinai J. Med. **78**, 280–290 (2011)

2. Tang, J., Rangayyan, R.M., Xu, J., Naqa, I.E., Yang, Y.: Computer-aided detection and diagnosis of breast cancer with mammography: recent advances. IEEE Trans. Inf. Technol. Biomed. **13**(2), 236–251 (2009). http://scholar.google.com/scholar?hl=en&btnG=Search&q=intitle:Computer-Aided+Detection+and+Diagnosis+of+Breast+Cancer+With+Mammography+:+Recent+Advances#0

3. Bezerra, L.A., Oliveira, M.M., Rolim, T.L., Conci, A., Santos, F.G.S., Lyra, P.R.M., Lima, R. C.F.: Estimation of breast tumor thermal properties using infrared images. Signal Process. **93** (10), 2851–2863 (2013). http://dx.doi.org/10.1016/j sigpro.2012.06.002

4. Kapoor, P., Prasad, S.V.A.V.: Image processing for early diagnosis of breast cancer using infrared images. In: The 2nd International Conference on Computer and Automation Engineering (ICCAE), vol. 3, pp. 564–566 (2010). http://ieeexplore.ieee.org/lpdocs/epic03/wrapper.htm?arnumber=5451827
5. Amri, A., Saidane, A., Pulko, S.: Thermal analysis of a three-dimensional breast model with embedded tumour using the transmission line matrix (TLM) method. Comput. Biol. Med. **41** (2), 76–86 (2011). http://dx.doi.org/10.1016/j.compbiomed.2010.12.002

6. Jiang, L., Zhan, W., Loew, M.H.: Modeling static and dynamic thermography of the human breast under elastic deformation. Phys. Med. Biol. **56**(1), 187–202 (2011)

7. Sanchez-Palencia, E.: Non-homogeneous Media and Vibration Theory, vol. 127. Springer, Berlin (1980). http://adsabs.harvard.edu/abs/1980LNP...127.....S

8. Bakhalov, N., Panasenko, G.: Homogenisation: Averaging Processes in Periodic Media. Kluwer Academic Publishers, Dordrecht (1989)

9. Penta, R., Ambrosi, D., Shipley, R.J.: Effective governing equations for poroelastic growing media. Q. J. Mech. Appl. Mech. **67**(1), 69–91 (2014)

10. Pennes, H.: Analysis of tissue and arterial blood temperatures in the resting human forearm. J. Appl. Physiol. **1**(2), 93–122 (1948)

11. Wahab, A.A., Salim, M.I.M., Ahamat, M.A., Manaf, N.A., Yunus, J., Lai, K.W.: Thermal distribution analysis of three-dimensional tumor-embedded breast models with different breast density compositions. Med. Biol. Eng. Comput. **1–11** (2015). http://dx.doi.org/10.1007/s11517-015-1403-7

12. Davit, Y., Bell, C.G., Byrne, H.M., Chapman, L.A.C., Kimpton, L.S., Lang, G.E., Leonard, K.H.L., Oliver, J.M., Pearson, N.C., Shipley, R.J., Waters, S.L., Whiteley, J.P., Wood, B.D., Quintard, M.: Homogenization via formal multiscale asymptotics and volume averaging: how do the two techniques compare? Adv. Water Resour. **62**, 178–206 (2013). http://dx.doi.org/10.1016/j.advwatres.2013.09.006

13. Muskhelishvili, N.I.: Some Basic Problems of the Mathematical Theory of Elasticity. Springer, Netherlands (1977)

14. Cionarescu, D., Donato, P.: An Introduction to Homogenization. Oxford University Press Inc., New York (1999)

15. Bardati, F., Iudicello, S.: Modeling the visibility of breast malignancy by a microwave radiometer. IEEE Trans. Biomed. Eng. **55**(1), 214–221 (2008)