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

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++.


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.

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].
ij @u e @x j n i ¼ ÀK g ij @u g @x j n i on @X 2 ; where K g ij ¼ k g d ij denotes glandular tissue thermal conductivity, q 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 2 K [11]. Besides, the rapidly oscillating coefficients K e ij ; x e b and q e m are defined as follows and q e m ðxÞ ¼ 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 X 1 with Lipschitz boundary @X 1 ¼ @X n 1 [ @X d 1 , where @X n 1 \ @X d 1 ¼ ;. On the other hand, the cancerous tissue will be characterized by two regions of dissimilar thermal properties: the tumoral area (X e t -fibers) and the glandular area (X e g -matrix). In this sense, the cancerous region will consist of a periodic microstructure associated with the open, bounded, and connected domain X 2 ¼ X e g [ X e t [ @X e g with Lipschitz boundary @X 2 ¼ @X e g and with X e g \ X e t ¼ ;. Then, the cancerous breast is represented by X ¼ X 2 [ X 1 (Fig. 1). Let e [ 0 be the size of the microstructure and y ¼ x=e 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 @Y t such that . It is also assumed that X e g is connected and that the inclusions do not intersect the boundary @X e g . In previous works, soft tissues assume to present this type of arrangement. In fact, Penta et al. [9] 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 @X n 1 and a prescribed temperature on @X d 1 . In the case of (2) we assume heat and temperature continuity on @X 2 . Moreover, continuity conditions for temperature and heat flow are imposed on C e (boundary between the glandular tissue X e g and the tumor inclusions X e t ), i.e., 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 e [problem (2)] is sought as a function of e for e ! 0, namely where the functions v p ðyÞ; u 2 ðx; yÞ, are Y-periodic in y. In particular, the vector function vðyÞ satisfies the unit cell problem À @ @y i K ij ðyÞ @v p ðyÞ and u 0 ðxÞ the homogenized problem solution where j Á j represents volume fraction. The effective constant coefficientsK ip are given byK where p ¼ 1; 2 and hÁi denotes volume average.

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 and v ðtÞ 2 ¼ Im where ðcÞ with c ¼ g; t denotes the constituent, the superscript o specifies that the sum is carried out over odd indices, the coefficients a p 0 ; a p k and c p l ðp ¼ 1; 2Þ are real and f is the zeta quasi periodic Weierstrass function. Now, using Laurent's expansion of f and the quasi-periodicity property of f and its derivatives v ðgÞ 1 ¼ Re where for p ¼ 1; 2 5 and S k are called the reticulate sums and are defined as S k ¼ P w2LÃ 1 w k ðk ! 3; k oddÞ with L Ã representing the lattice excluding the number w ¼ 0 and w ¼ mw 1 þ nw 2 where m; n 2 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 C; z ¼ Re ih where R is the circumference radius give that coefficients a p k 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) Using the form of K ij , Green's theorem, the double periodicity of v p and formulas (8)-(10), thenK In fact, if k g ¼ k t . Then,K ¼K 11 ¼K 22 .

Numerical Solution and Analysis of Results
This section is devoted to find the temperature distribution of problems (1) and (6) where we define as With this aim, we follow the following procedure.

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 1 d ðX 1 Þ ¼ fu 2 H 1 ðX 1 Þ s:t: cðuÞ ¼ 0 on @X d 1 g. Using the trace theorem for u c 2 H 1=2 ðX 1 Þ, there exists a continuous linear operator R 0 : In this way, the equivalent variational formulation of problem (1) is where aðũ; vÞ ¼ hðu e À R 0 u c ÞvdS: 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ðũ; vÞ is continuous In this sense, observe that K g 2 Mða; b; X 1 Þ and by Cauchy-Schwartz Furthermore, by the Poincaré-Friedrichs II theorem On the other hand, by the trace theorem Therefore, (ii) The bilinear form aðũ; vÞ is H 1 d -elliptic Let u 2 H 1 d ðX 1 Þ, (iii) The linear form LðvÞ is continuous in But, by Cauchy-Schwarz inequality, the Poincaré-Friedrichs II theorem and the fact that R 0 u c 2 H 1 ðX 1 Þ and r x ðR 0 u c Þ 2 L 2 ðX 1 Þ ð Þ n jLðvÞj C 10 kvk H 1 d ðX 1 Þ ; where Then, 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 r x u Á n þ u=e ¼ u c =e on @X d 1 and solves the problem with a very small value of e. In particular, we approximate the involved functions by piecewise linear continuous finite elements.

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 1 0 ðX 2 Þ ¼ fu 2 H 1 ðX 2 Þ s:t: cðuÞ ¼ 0 on @X 2 g. Using the trace theorem for u g 2 H 1=2 ðX 2 Þ there exists a continuous linear operator R 0 : H Then, on @X 2 cðũ 0 Þ ¼ cðu 0 Þ À cðR 0 u g Þ ¼ u g À u g ¼ 0: In this way, the equivalent variational formulation of problem (6) is where aðũ 0 ; vÞ ¼ hg e iðR 0 u g Þvdx: In particular,K 2 Mða; b; X 2 Þ see Cionarescu and Donato [14], and hf e i 2 L 2 ðX 2 Þ and hg e i [ 0. The existence and uniqueness of solutionũ 0 can be proved through the Lax-Milgram theorem. Then, we must show that: (i) The bilinear form aðũ 0 ; vÞ is continuous In this sense, observe that using the fact thatK 2 Mða; b; X 2 Þ and by Cauchy-Schwartz inequality jaðũ 0 ; vÞj bkr xũ 0 k L 2 ðX 2 Þ kr x vk L 2 ðX 2 Þ þ hg e ikũ 0 k L 2 ðX 2 Þ kvk L 2 ðX 2 Þ : Now, from remark 3.37 p. 32 by Cionarescu and Donato [14], for u 0 ; v 2 H 1 0 ðX 2 Þ, Furthermore, by the Poincaré-Friedrichs I theorem, (iii) The linear form LðvÞ is continuous in jhg e iðR 0 u g Þvjdx: But, using the Cauchy-Schwarz inequality, the Poincaré-Friedrichs I theorem and the fact that R 0 u g 2 H 1 ðX 2 Þ and r x ðR 0 u g Þ 2 L 2 ðX 2 Þ ð Þ n jLðvÞj C 6 kvk H 1 0 ðX 2 Þ ; where C 6 ¼ C 4 khf e ik L 2 ðX 2 Þ þ hg e iC 5 kR 0 u g k L 2 ðX 2 Þ þ bkr x ðR 0 u g Þk L 2 ðX 2 Þ : Thus, (i)-(iii) proves the existence and uniqueness of solutionũ 0 by using the Lax-Milgram theorem. Now, it must be shown that the map hf e i 2 L 2 ðX 2 Þ ! u 2 H 1 0 ðX 2 Þ is continuous in order to prove the regularity of the weak solution. In fact, from the H 1 0 -ellipticity of the bilinear form jaðu; uÞj ! C 3 kuk 2 Then, i.e., 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.

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 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 t m ¼ C=ð468:6 lnð100DÞ þ 50Þ, where C ¼ 3:27 Â 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 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 jY t j ¼ 0:7 is considered so that healthy breast tissue volume fraction is jY g j ¼ 0:3 in the tumorous region. As observed, if the breast dimensions are bigger, the maximum temperature is higher,     [15] 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.

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 then 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.