The heat transfer analysis of nanoparticle heat source in alanine tissue by molecular dynamics

(1)

The heat transfer analysis of nanoparticle heat source in

alanine tissue by molecular dynamics

David T.W. Lin

a

, Ching-yu Yang

b,∗

a_{Department of Information Management, Hsing Kuo University of Management, No. 89, Yuying St., Tainan City 709, Taiwan, ROC} b_{Department of Mold and Die Engineering, National Kaohsiung University of Applied Sciences,}

415 Chien Kung Road, Kaohsiung City 807, Taiwan, ROC

Received 1 May 2005; received in revised form 27 June 2005; accepted 27 June 2005

Abstract

The purpose of this study is to simulate the heat transfer problem when the 3-D Alanine tissue is heated by the gold nanoparticle in the field of molecular dynamics. In this paper, the Alanine molecule is adopted and its parameters are available in the GROMACS protein data bank. A computing algorithm is developed to evaluate the heat transfer phenomena in the nano-scale biological system based on the molecular dynamics and the protein data bank. The value of the thermal conductivity of Alanine is calculated from the autocorrelation function of the Green–Kubo formalism and this result has a roughly approximation with the bulk thermal conductivity reported by experimental data[1]. Two kinds of problems are investigated in the paper. One is the Alanine tissue heated by the constant heat source and the other is by the time-varying heat source. The numerical results show that a temperature jump exists around the source and the temperature profiles drop to the environmental temperature within a very short distance. It concludes that only a small region around the nano-scale heat source is affected by the heated process. Therefore, the results of the nanoparticle-heated method could be applied to the clinical therapy of tumor, and the normal cells are destroyed only within a smaller region than those of chemotherapy or surgery.

Keywords: Nanoparticle; Molecular Dynamics; Alanine; Protein data bank; Tumour

1. Introduction

During the past decade, a new thermal therapy was devel-oped to treat the pathological cells directly through the nanoparticles such as the magnetic fluid hyperthermia (MFH)

[2] and the therapy of light to heat [3]. Additionally, the nanoparticles coated with the suitable antibodies are used as labels to detect the corresponding tumor cell and to adsorb the surface of the tumor cell. Many nanoparticles will respond to an external applied field, i.e., magnetic field or focused light and so on. The nanoparticles convert the absorbed energy of the external field to heat and destroy the adsorbed cell. According to the mechanisms of the hyperthermia and

ther-∗_{Corresponding author. Tel.: +1 886 7 3814526 5413;} fax: +1 886 7 383 5015.

E-mail address: [email protected] (C.-y. Yang).

moablation, the death of the cell is expected[4]. Therefore, the tumor cells are destroyed under the controllable power of the external field. Consequently, it is important to understand the heat transfer phenomena included the temperature pro-file and the thermal conductivity in the nano-scale bio-tissue. Nevertheless, the classical continue theory is not befitted to the nano-scale phenomena because the particle’s motion becomes the major effect in the nano-scale system. As well, some macro-scale phenomena are not existed in the nano-scale such as the non-slip condition and the local heat equi-librium. Some special phenomena in the nano-scale problem are the thermal jump and no natural convection because the size effect.

Redondo and Lesar [5] studied the bio-material in the

atomistic level. According to their report, the molecular dynamics and the Monte Carlo method are the major method-ology to evaluate the properties of the molecules and

(2)

blies of the different molecules. That is the reason that the molecular dynamics is chosen to treat the nano-heat transfer problem in this paper. In addition, Hilger et al.[4]presented the elimination of small tumours in the breast using MFH fur-ther. The approach reviewed the biophysical basis, the target temperature control and the degree of energy tuning. Flen-niken et al.[6]demonstrated that the small heat shock protein cage is a new nano-scale platform whose exterior and interior surfaces are amenable to both genetic and chemical modifi-cation.

The bio-tissue material is association of cells that is con-structed by protein and aminoacid. Therefore, this paper studies the heat transfer phenomena of the aminoacid and lead to understand the nano-scale bio-tissue thermal problem. The elevation of the protein’s force field is the key process in the molecular dynamics and the force field is very complex. Many previous studies have been evaluated the properties of the protein by the protein data bank in the molecular dynam-ics. Bastawissy et al.[7]proposed the research included the protein material at normal and elevated temperature in the molecular dynamics. Suenaga[8]investigated a small-sized protein folding with implicit solvent. Hornak and Simmerling

[9] developed the softcore potential functions to overcome steric barriers. Yamashita et al.[10]investigated the compar-ison of sampling efficiency between the molecular dynamics and the Monte Carlo method in protein. Komeiji et al.[11]

studied the protein simulation by parallel molecular

dynam-ics. Zal and Gascoigne[12]used the live FRET imaging to

reveal early protein-protein interactions.

Many researchers investigated the dynamic transport properties of the different materials in order to build the roadmap of the nano-scale thermal sciences. The

Green–Kubo formalism [13]is often used to calculate the

dynamic transport properties such as the thermal conductiv-ity, thermal diffusivity and so on. Andrade and Stassen[14]

calculated the thermal conductivity from the Green–Kubo time correlation function. Fren´andez et al.[15]investigated the shear and bulk viscosity and thermal conductivity of sim-ple real fluids. The results have good agreement with the experimental data except the bulk viscosity. According to the results of McGaughey and Kaviany, the thermal conductivity is temperature dependent or independent that is affected by how long the time scale[16]. For the comparison of the results of the thermal conductivity of the protein, some previous researches in the protein bulk thermal conductivity have been reviewed. Karen et al.[1]presented the optical and thermal characterization of albumin protein solders and the thermal conductivity of protein solder was from 0.36 to 0.4 W/m K. The materials enriched proteins are studied by Krokida et al.[17]and the thermal conductivity of these samples were 0.3–0.6 W/m K.

The purpose of this research is to investigate the basic phenomena of the nano-scale heat transfer problem in order to propose a bio-tissue heat transfer model. This study treats the characteristics of nano-scale heat transfer in a nano-cavity with the consideration of the realize interactions of

macro-molecules. In the calculated process, the combination of the protein data bank and the methodology of molecular dynam-ics is used to develop a new algorithm to evaluate the dynamic transport properties in the nano-sized biological material. The results can be applied in the clinical recognition of tumour therapy and the non-invasive photo-thermal tumours ablation diagnosis.

2. Mathematical model

In this paper, the heat transfer problems in Alanine are examined numerically by the molecular dynamics with the leap-frog method and the GROMACS protein data bank. The purpose of this study is to discuss the evolution of the temper-ature in the thermostate, the transient models of the Alanine and the thermal conductivity of the Alanine when the system is heated by the gold nanoparticle.

The interaction phenomenon between the molecules is an important issue to investigate the physical objects’ state, i.e. solid, liquid and gas. Furthermore, the understanding of the molecular interaction is the way to study the dynamic prop-erties of the physical objects. In the past, the behavior of the molecular particle is difficult to trace due to the lim-its of the computer capability when Lagrangian’s theorem is implemented. However, it is possible to trace the interac-tions between molecules today when the molecular dynamics theory is applied. Molecular dynamics is used to solve the solution of the classical equations of motion (Newton’s equa-tions) when the interaction of a set of molecules is considered.

fi= −∇rφ(rij) (1)

whereφ is the potential between the molecules i and j, fiis the force acted on molecule i.

A leap frog method (i.e., modified Verlet algorithm) is modified from the basic Verlet scheme with a half-step scheme.

The algorithm is formulated as follows:

¯ r(t + t) = ¯r(t) + t ¯v t −1 2t (2) ¯ v t +1 2t = ¯v t −1 2t + t ¯a(t) (3)

where ¯r,¯v, and ¯a(t) are the position, velocity, and acceleration of each molecule at each time stept.

According to Eqs.(1)–(3), the position, velocity and accel-eration of each molecule at each time step can be calculated. In this research, the Andersen thermobath method is used to simulate the interaction between the Alanine

macro-molecule and the gold nanoparticle [18]. Moreover, the

molecular velocity is corrected at each time step when the system is maintained in a particular temperature and the

(3)

cor-rection of molecular velocity is listed as below: T_a∗= 1 3N i v∗_i2 , v∗,new_i = v∗_i T∗ s T∗ a (4)

where N is numbers of water molecules in the simulated sys-tem,T_a∗ is the instantaneous dimensionless temperature of the simulated system after all colliding processes at each time step, T_s∗ is the initial dimensionless temperature of the simulated system, v∗_i is the dimensionless velocity of the ith molecule after a colliding process at each time step, andv∗,new_i is the dimensionless corrected molecular velocity of the ith molecule, based on the statistical thermodynam-ics, it is known that the initial velocities of the molecules in the equilibrium system are distributed according to the Maxwell–Boltzmann velocity distribution while the temper-ature of the isolated system is fixed.

In the transient case, the net energy follows the ther-mobath process in the interaction between nanoparticle and the Alanine macromolecules. Therefore, the temperature of the nanoparticle changes with the variation of the Alanine’s kinetic energy. The boundary of the cavity follows the ther-mobath condition and the temperature of the therther-mobath is fixed at the environmental temperature. Therefore, the tem-perature of the nanoparticle decays as the runs. The formu-lations describe as below:

mpcp(T_t_i+1_,p− T_t_i_,p)= −m 2 j (V_t2 i+1,j− V 2 ti,j) (5)

whereV_t_i_,jis the velocity of Alanine j at time ti.Vti,pis the

temperature of the gold nanoparticle at time ti. m and mpare

the mass of Alanine and the gold nanoparticle. As well, cpis

the specific heat of the gold nanoparticle.

In addition, the net momentum of the system must be van-ished in order to guarantee that the system will not move due to an external force. The conservation conditions for the momentum are described as below.

v∗,new_i = v∗ i − 1 N N i=1 p∗ i (6)

wherep∗_i is the dimensionless momentum of the ith molecule. In the Green–Kubo method, the thermal conductivity k is related to the time integral over the heat current autocorrela-tion funcautocorrela-tion. The funcautocorrela-tion of k is shown as follows[14]:

k = 1

kBVT2

_∞

0 Jx

(0)J_x(t) dt (7)

where kB, V, and T are Boltzmann’s constant, system’s

vol-ume and temperature of system. The heat current J is:

J = 1 2  N i=1 mv2 iv¯i+ N i=1 N j =i φ(rij)¯v_i− w(r_ij)r¯ij rijv¯i ¯ rij rij   (8)

whereφ(rij) is the potential between molecule i and j. Addi-tionally,w(r_ij) is the pair virial function and defines as:

w(rij)= r_ijdφ(rij)

dr_ij (9)

3. Simulation details

The temperature of the heat source (nanoparticle) is mod-eled as time-independent (thermostate) and time-dependent (transient), separately. The initial temperature of the heat source (Tini) with the thermostate cases are 473, 423, 373 and

323 K, respectively. The Tiniis assumed as 373 K in the

tran-sient case and the total energy of the heat source is decreased while the interaction between the Alanine macromolecule and the heat source is active. Therefore, the temperature of the heat source gradually decreases to the environmental temperature when the computing time increases. The num-ber of the Alanine is 345. The cross-sections of this system

in x–y plane are shown in Fig. 1. The size of the system

(X× Y × Z) are 700 nm × 700 nm × 700 nm, i.e., aspect ratio is unity, and its coordinates are−X/2 to X/2, −Y/2 to Y/2 and

−Z/2 to Z/2. It implies that the global density (ρ) is about

1000(N/␮m3) in all cases. The diameter of heat source is 100 nm and the source locates at the center of this system. The chemical formula of Alanine used in this research is NH C2H4CO. The topology of the Alanine is constructed from

the building block of the GROMACS protein data bank. The GROMACS protein data bank is used to model the compli-cated interaction between two Alanine macromolecules[19]. The position, velocity and acceleration of each molecule are exported from the GROMACS and these converted to the new initial data for succeed simulation by a half-step modified Verlet algorithm[20,21]. In the molecular dynamics, the trun-cation distance of potential must be considered. According to the research of Allen and Tildesley[22], the dimension-less cutoff distance is r∗_c = l/2 (i.e., l is the dimensionless length of the channel) when the complicated potential is used.

The simulations are started from the equilibrium state of an initially crystalline fluid. The velocity distribution is satisfied with the algorithm of Maxwell velocity distribution and a similar period of equilibration with the switched-on external field. In this work, the total time steps are 2.5× 107with a time step (t = 0.1 fs) except the equilibrated process. Ten thousand independent autocorrelation functions are used to calculate the thermal conductivity. According to Schlick’s study [23], the stable criterion for of the leap-frog method ist = Γ /π for a harmonic oscillator of period, Γ . The time intervalt is 0.1 fs that is much less than 6.4 fs for H O H bend and 3.1 fs for O H stretch. All boundary conditions are periodic boundary condition (PBC) in the process. The initial temperature is T = 309 K and the pressure of the system and the environment is P = 1 atm. The parameters of the simulated model are listed inTable 1.

(4)

Fig. 1. The cross-section of simulation geometry on x–y plane.

4. Computational algorithm

The proposed iterative method is summarized as follows and inFig. 2:

Step 1: The GROMACS protein data bank is used to build the initial model and the sample molecule is Ala-nine macromolecule. A regular triangular grid in the plane is constructed by the initial position of each molecule. During a period of time-steps, the system relaxes from the initial conditions and it approaches to the

equilib-rium state. The equilibequilib-rium process is shown in steps 2–5.

4.1. Equilibrium process

Step 2: Through the GROMACS protein databank, the acceleration ¯a_iis solved to represent the interaction between two Alanine molecules.

Step 3: The content of the GROMACS output file includes the position, the velocity and the acceleration of the molecule at each time step. This file is treated as a new

Table 1

The parameters of the simulated model

Temperature of the enviroment 309 K

Initial temperature of the system 309 K

Pressure 1 atm

Size of the system 700 nm× 700 nm × 700 nm

Diameter of the nanoparticle D = 100 nm

The sample of molecule Standard alanine macromolecule

Protein data bank GROMACS

Simulated time 2.5× l07_runs _{Time step}_t _0.1fs

Simulated problem: (1) thermo-state problem The initial temp. of the nanoparticle 473, 423, 373, and 323 K; (2) transient problem. The initial temperture of the nanoparticle 373 K.

(5)

(6)

initial input file for the program to calculate the intermedi-ated position, the intermediintermedi-ated velocity of each molecule by the leap-frog method.

Step 4: The velocity of each molecule is corrected through Eq.(4)at every time step to maintain the temperature of the system.

Step 5: If the velocity distribution is satisfied with the algo-rithm of Maxwell velocity distribution then the process is terminated, otherwise the computation is returned to step 2. 4.2. Start simulation

Step 6: Repeat the process of step 2 to obtain the accelera-tion.

Step 7: Repeat the process of step 3 to obtain the position and velocity of each molecule.

Step 8: The thermobath is considered to correct the velocity of each molecule and to build the temperature profile. Step 9: In order to maintain the temperature and the momentum conservation of the system, the velocity of each molecule is corrected through Eqs.(4)and(6)at each time step.

Step 10: Terminate the process if the computed time approaches, otherwise the computation is returned to step 6.

5. Results and discussions

In this paper, two cases are investigated. In the first case, the thermostate case is considered. The Alanine molecules are heated by the fixed heat source in the cavity. The Tini

is assumed as 473, 423, 373 and 323 K, respectively. The temperature profiles at x = 0–350 nm, y = 0 nm, z = 0 nm are presented inFig. 3. A close examination ofFig. 3discloses that the temperature profiles have the temperature jump close to the heat source. The temperature of heat source is 473 K and the width of the temperature jump region is about 20 nm (x = 100–120 nm). The temperature is from 473 K down to

Fig. 3. The temperature profiles of the thermostate simulation.

453 K in the temperature jump region. Behind this region the classical linear temperature profile is observed and it obtained by the Fourier conduction equation when the distance is over jump region (i.e., x = 100–350 nm). It is clear that the temper-ature profiles at Tini= 423, 373 and 323 K are similar to that

of Tini= 473 K. In general, the temperature profile is divided

into two regions. One is the boundary jump and the other is the linear drop. Nevertheless, the jump boundary in this research is different from the classical temperature profiles. The phe-nomena of the boundary jump is dominated in the nano-scale problem but it is minor in the macro-scale problem. This rea-son is that the nano-scope boundary is not able to describe a proper local equilibrium due to the complicated interaction of the approaching molecules and the leaving molecules at the boundary. The similar local effect (i.e., slip condition) is also presented by the previous research[24]. Therefore, it is imperative to derive a proper boundary for the nano-scope problem. The maximum diameter of the region which is affected by the heat source is about 550 nm at Tini= 373 K

(seeFig. 3). The number of the Alanine molecules are affected by the heat source is about 60 and the average diameter of the macromolecule is assumed as 200 nm. Therefore, the normal cells destroyed through the light-to-heat therapy are smaller than the chemotherapy or surgery.

The temperature profile with a transient heat source is pre-sented inFig. 4. The temperature of Tiniis started at 373 K.

The internal energy of the heat source is dissipated to envi-ronment when the time increases gradually. Both the heat source temperature and the system temperature are changed to the environmental temperature when the action is finished. The evolution of the temperature profile is shown inFig. 4

when x = 100 nm, y = 0 nm, z = 0 nm. The result shows that the system temperature increases to the temperature peak first. When the maximum temperature reaches, the system temper-ature decreases to the environmental tempertemper-ature gradually. Both the nano-scale and macro-scale temperature profiles have the same trend in the transient problem. Moreover, the

Fig. 4. The temperature profile of the transient simulation at x = 100 nm, y = 0 nm.

(7)

Fig. 5. The thermal conductivity profile of Alanine in the thermostate case.

system temperature decreases to the environmental tempera-ture quickly.

Finally, the thermodynamic property of nanoscale-Alanine is calculated in this paper,Fig. 5shows the behavior of the thermal conductivity for the Alanine in the thermostate case. The thermal conductivity rises to the stable value after 135ps. The value of thermal conductivity is approximated to 0.42 W/m K either the temperature is high or low. In addi-tion, the numerical result shows that the estimation of the thermal conductivity has an agreement with that of the bulk protein property. Consequently, the thermal conductivity has the same order in the in the macro- or nano-scale bio-heat transfer problem.

6. Conclusion

The phenomena of the heated Alanine tissue by gold nanoparticle are investigated in this paper. The GROMACS protein data bank is used to predict the Alanine temperature profiles and the thermal conductivity is solved based on the leap-frog method of the molecular dynamics. In this research, the temperature profiles have the significant jump near the heat source that is deviated from the result calculated from the classical heat conduction theory. As well, a tiny heated region is affected by the heat source. It is also shown that the system temperature recovers to the environmental tempera-ture within a very short time. In concluded that the maximum system temperature and the heated region are controllable as long as the temperature of the heat source is controllable. In other words, a good control of the nanoparticle energy level leads to a precision bio-tissue heating in the clinical ther-apy. Moreover, the value of the thermal conductivity is about

0.42 W/m K. The thermal conductivity of the Alanine macro-molecule is firstly derived from the Green–Kubo formalism in this paper. Moreover, this paper constructs a bio-tissue heat transfer model and it expands the application of thermal science to the bio-medical field.

Acknowledgement

The computer time on the IBM SP2 supercomputers of this work by the National Center for High Performance Comput-ing of Taiwan, ROC is gratefully acknowledged.

References

[1] M.M. Karen, S.S. Brian, C.B. Naresh, G.D. Mathieu, J.W. Ashley, M.D. Judith, Appl. Optic. 38 (31) (1999) 6661–6672.

[2] A. Jordan, R. Scholz, P. Wust, H. Fahling, R. Felix, J. Magn. Magn. Mater. 201 (1999) 413–419.

[3] D.P. O’Neal, L.R. Hirsch, N.J. Halas, J.D. Payne, J.L. West, Cancer Lett. 109 (2004) 176–181.

[4] I. Hilger, R. Hergt, W.A. Kaiser, IEEE Proc. Nanobiotechnol. 152 (1) (2005) 33–39.

[5] A. Redondo, R. Lesar, Annu. Rev. Mater. Res. 34 (2004) 279–314. [6] M.L. Flenniken, D.A. Willits, S. Brumfield, M.J. Young, T. Douglas,

Nano Lett. 3 (11) (2003) 1573–1576.

[7] E. El-Bastawissy, M.H. Knaggs, I.H. Gilbert, J. Mol. Graph. Model. 20 (2001) 145–154.

[8] A. Suenaga, J. Mol. Struct. (Theochem.) 634 (2003) 235–241. [9] V. Hornak, C. Simmerling, J. Mol. Graph. Model. 22 (2004)

405–413.

[10] H. Yamashita, S. Endo, H. Wako, A. Kidera, Chem. Phys. Lett. 342 (2001) 382–386.

[11] Y. Komeiji, M. Haraguchi, U. Nagashima, Parallel Comput. 27 (2001) 977–987.

[12] T. Zal, N.R.J. Gascoigne, Curr. Opin. Immunol. 16 (2004) 418– 427.

[13] D.A. McQuarrie, Statistical Mechanics, Harper Collins, New York, 1976.

[14] J. de Andrade, H. Stassen, J. Mol. Liquids 110 (2004) 169–176. [15] G.A. Fern´andez, J. Vrabec, H. Hasse, Fluid Phase Equilib. 221

(2004) 157–163.

[16] A.J.H. McGaughey, M. Kaviany, Int. J. Heat Mass Transfer 47 (2004) 1783–1798.

[17] M.K. Krokida, P.A. Michailidis, Z.B. Maroulis, G.D. Saravacos, Int. J. Food Propert. 5 (1) (2002) 63–111.

[18] H.C. Andersen, J. Chem. Phys. 72 (1980) 2384–2393.

[19] H.J.C. Berendsen, D. van der Spoel, R. van Drunen, Comput. Phys. Commun. 91 (1995) 43–56.

[20] R.W. Hockney, Methods Comput Phys. 9 (1970) 136–211. [21] D. Potter, Computational Physics, Wiley, New York, 1972. [22] M.P. Allen, D.J. Tildesley, Computer Simulation Liquids, Clarendon

Press, Oxford, 1987.

[23] T. Schlick, E. Barth, M. Mandziuk, Annu. Rev. Biophys. Struct. 26 (1997) 181–222.