Journal of Colloid and Interface Science 437 (2015) 187–196

Contents lists available at ScienceDirect

Journal of Colloid and Interface Science www.elsevier.com/locate/jcis

Molecular dynamics simulations of proton transverse relaxation times in suspensions of magnetic nanoparticles Tomasz Panczyk a,⇑, Lukasz Konczak a, Szczepan Zapotoczny b, Pawel Szabelski c, Maria Nowakowska b a

Institute of Catalysis and Surface Chemistry, Polish Academy of Sciences, ul. Niezapominajek 8, 30239 Cracow, Poland Faculty of Chemistry, Jagiellonian University, ul. Ingardena 3, 30060 Cracow, Poland c Department of Chemistry, Maria Curie-Sklodowska University, pl. M. Curie-Sklodowskiej 3, 20031 Lublin, Poland b

a r t i c l e

i n f o

Article history: Received 3 May 2014 Accepted 29 August 2014 Available online 19 September 2014 Keywords: Magnetic nanoparticles T2 relaxation Molecular dynamics Intermolecular interactions

a b s t r a c t In this work we have analyzed the influence of various factors on the transverse relaxation times T2 of water protons in suspension of magnetic nanoparticles. For that purpose we developed a full molecular dynamics force field which includes the effects of dispersion interactions between magnetic nanoparticles and water molecules, electrostatic interactions between charged nanoparticles and magnetic dipole–dipole and dipole–external field interactions. We also accounted for the magnetization reversal within the nanoparticles body frames due to finite magnetic anisotropy barriers. The force field together with the Langevin dynamics imposed on water molecules and the nanoparticles allowed us to monitor the dephasing of water protons in real time. Thus, we were able to determine the T2 relaxation times including the effects of the adsorption of water on the nanoparticles’ surfaces, thermal fluctuations of the orientation of nanoparticles’ magnetizations as well as the effects of the core–shell architecture of nanoparticles and their agglomeration into clusters. We found that there exists an optimal cluster size for which T2 is minimized and that the retardation of water molecules motion, due to adsorption on the nanoparticles surfaces, has some effect in the measured T2 times. The typical strengths of the external magnetic fields in MRI are enough to keep the magnetizations fixed along the field direction, however, in the case of low magnetic fields, we observed significant enhancement of T2 due to thermal fluctuations of the orientations of magnetizations. Ó 2014 Elsevier Inc. All rights reserved.

1. Introduction Superparamagnetic colloids improve the contrast in Magnetic Resonance Imaging (MRI). Thanks to their great magnetization, the nanoparticles induce great variations in the surrounding water-proton relaxation times. Among many possible materials for fabrication of magnetic nanoparticles [1] iron oxide is the most promising since it is nontoxic, and fabrication of nanoparticles does not require hazardous organic solvents and extreme conditions of pressure or temperature [2,3]. Iron oxide nanoparticles (NPs) predominantly shorten T2 transverse relaxation time and thus provide negative contrast in T2 weighted images. Moreover, iron oxide NPs can easily be coated by protective layers enhancing their colloidal stability and biocompatibility. The T2 relaxation time depends on many parameters of the NPs, therefore their function as contrast agents is also a function of ⇑ Corresponding author. Fax: +48 81 537 5685. E-mail address: [email protected] (T. Panczyk). http://dx.doi.org/10.1016/j.jcis.2014.08.066 0021-9797/Ó 2014 Elsevier Inc. All rights reserved.

those parameters. Thus, a careful analysis of those parameters allows for recognizing optimal ranges of their values and designing synthesis procedures leading to possibly best outcomes. The decrease in T2 has been explained theoretically for single spherical NPs by the inner- and outer-sphere theory [4], then refined with chemical exchange [4,5] and partial refocusing models [6]. When theoretical models became invalid, particularly for analysis of NPs agglomerates, computational methods were employed. Those methods were based on a simple random walk scheme for water-protons motion while NPs were normally static. Also, no intermolecular interactions (except of hard-core repulsion) were taken into account [6–10]. Though simple, the random walk methods allowed drawing very useful conclusions concerning the NPs sizes, agglomeration, cluster sizes and structure in the modification of the transverse relaxation times. The aim of this work is the analysis of various parameters of magnetic NPs in the modification of T2 relaxation time but we focus on possibly a complete description of intermolecular interactions and the dynamics of motion of both water-protons and magnetic nanoparticles. For that purpose we designed a full molecular

188

T. Panczyk et al. / Journal of Colloid and Interface Science 437 (2015) 187–196

dynamics model of the NPs suspension in aqueous media accounting for dispersion forces, electrostatic interactions and magnetic interactions including magnetization reversal within the nanoparticle volume. All those factors are carefully studied and conclusions concerning their meaning are provided. Our calculations cover the range of parameters representative of chitosan covered iron oxide nanoparticles. Those NPs have recently been synthesized and characterized in our group [3]. It was found that they form highly stable and biocompatible suspensions; moreover, they reveal very high relaxivity values. Thus, another aim of this work is devoted to theoretical understanding of the observed properties of chitosan covered nanoparticles and finding possible routes for optimization of their properties. 2. Methods All calculations were performed using the Large-scale Atomic/ Molecular Massively Parallel Simulator (lammps) code [11] with several extra classes written from scratch and working on magnetic torques. The force field consisted of several interactions types which are shortly discussed in subsequent paragraphs. 2.1. Magnetic Interactions A superparamagnetic nanoparticle is modeled as a sphere of radius r = rM + d, that is, it is composed of a magnetic core of radius rM and nonmagnetic shell of thickness d. The magnetic core is characterized by some value of the saturation magnetization Ms, thus the particle magnetization

~ ¼ M s V½mx ; my ; mz  m

ð1Þ

where mi are unit vectors of magnetization and V ¼ ð4=3Þpr 3M . Each NP has also an easy axis of magnetization ~ e which defines an equilibrium direction of magnetization. In the absence of external magnetic field the magnetization aligns with the direction of the easy axis and any displacement of the magnetization is accompanied by the anisotropy energy increase 2

Ea ¼ K a sin h

ð2Þ

~ ~ e and Ka is the magnetic anisotropy constant. where cosh ¼ m Thus, any displacement of magnetization orientation induces the anisotropy torque ~ sa acting on the magnetization vector

~ Ea ~ ~ r sa ¼ m

ð3Þ

The same torque but of opposite direction acts on the easy axis vector. The presence of other magnetic NPs and external magnetic field, B, induce magnetic torques due to dipole–dipole interactions ~ sd and dipolefield interactions ~ sB . The dipole–dipole interactions can be adequately described by the Dormann–Bessais–Fiorani model [12] which yields the following formula for dipole–dipole torque

"

~i  ~ ~j  ~ ~i  m ~j dÞðm dÞ m l 3ðm ~ sd ¼ 0  5 3 4p d d

# ð4Þ

where l0 is the magnetic permeability of free space and d is the distance between two dipoles. The dipole–field interaction leads to ~ sB torque acting on the dipole

! ~ ~ B sB ¼ m

ð5Þ

where B is the external magnetic flux density. The sum of the above contributions drives the Néel rotation (rotation of the magnetic moment between easy axes) of the NP magnetization, however, due to the coupling in Eq. (3) the Brown-

ian rotation (rotation of the NP as a whole within its body frame) occurs as well, that is

~ NP dx ¼ I1~ sa dt

ð6Þ

where I is the NP inertia. As seen, the Eq. (6) defines the angular acceleration of the Brownian rotation. It is driven by the magnetic anisotropy torque only and it is calculated from Eq. (3). The negative sign means that the direction of the torque acting on the easy axis is opposite to the torque acting on the magnetization vector. The instantaneous orientation of the easy axis can be found by integration of the following equation

d~ e ~ NP  ~ ¼x e dt

ð7Þ

In the case of fine magnetic nanoparticles revealing moderate anisotropy constants the inertialess Néel rotation is much faster than the Brownian rotation. Thus, the full description of the dynamics of magnetization reversal according to the Landau–Lifshitz–Gilbert equation [13,14] is not necessary and we may rely only on its stationary solution. That solution has been found by analyzing the Fokker–Planck equation governing the time evolution of the nonequilibrium probability distribution of magnetic moments orientations associated with the stochastic Landau–Lifshitz–Gilbert equation [13]. Thus, on introducing the appropriate Néel time t0, which is the characteristic time of free diffusion in the absence of potential, the stationary solution of the Fokker– Planck equation has the Boltzmann distribution when

t0 ¼

~j 1 jm k 2ckT

ð8Þ

where c is a gyromagnetic ratio and k is a dimensionless damping coefficient that measures the magnitude of the relaxation (damping) term relative to the gyromagnetic term in the dynamical equation. We can further assume that magnetization reversal proceeds according to a coherent rotation mechanism, that is, the magnitude of magnetization does not change during rotation of the moment and thermal agitation helps crossing the energy barrier associated with magnetic anisotropy (Néel–Brown model). That mechanism has been experimentally confirmed for the case of cobalt and some other nanoparticles [15]. So, the magnetization reversal is in our case a simple free rotational diffusion process with an activation barrier. Thus, a mean time spent at a given orientation tN follows the exponential law,

tN ¼ t0 expðE=kTÞ

ð9Þ

where E is the net activation barrier, 2

!

~ B: E ¼ K a V sin h  m

ð10Þ

Eq. (10) shows that the barrier for rotational diffusion is controlled mainly by the magnetic anisotropy energy and the energy of magnetic dipole in the external field. Formally, it also depends on the dipole–dipole energy but that component is negligible in our case. The dynamics of inertialess magnetization reversal of the magnetic moments is therefore described by the overdamped Langevin dynamics. According to that scheme the angular acceleration is zero because all forces (torques) and friction terms are compensated by stochastic terms and the motion becomes inertialess. It physically means that the Néel rotation becomes a simple rotational diffusion process. Under the above assumptions the angular ~ m reads, velocity of magnetization x

kc ~m ¼ ~ NP þ ~ x sm þ x Ms V

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2kckT ! R M s V Dt

ð11Þ

T. Panczyk et al. / Journal of Colloid and Interface Science 437 (2015) 187–196

189

where ~ sm ¼ ~ sa þ ~ sd þ ~ sB is the net magnetic torque acting on the ! ~ , Dt is integration timestep and R are uniformly distribmoment m uted random numbers with zero mean and a unit standard deviation. The Brownian and Néel rotations are, as seen, coupled by adding the angular velocity of the NP to the rhs of Eq. (11). The instantaneous orientation of magnetization is thus given by

shell thickness. However, it can be shown that for shell thicknesses larger than 5–10 Å interactions between cores are negligible compared to the interaction between shells. Therefore, such core–shell NPs behave like they would be composed entirely of the shell material [19,18].

~ dm ~m  m ~ ¼x dt

2.3. Electrostatic interactions

ð12Þ

The above set of equations allows us to model the dynamics of the magnetization reversal of NPs as a function of the anisotropy barrier and applied field. Though based on the simplest approaches to the treatment of the Néel rotation it preserves the most important properties of the coherent rotation model. 2.2. Dispersion interactions Magnetic NPs in solution form suspension of colloid nanoparticles. Interaction between individual NPs is a key factor which controls stability of the suspension. Because NPs are extended bodies the interparticle interactions can be adequately described using the Hamaker theory [16]. Thus, upon integration of contributions from individual atoms creating the NPs the effective pairwise interaction energy between two spherical particles of radius a is given by

U CC

" # 2 Acsc 2a2 a2 d  4a2 Acsc þ ¼ þ 2 2 þ ln 2 6 d2  4a2 37800 d d " # 2 2 6 2 r d  14ad þ 54a d þ 14ad þ 54a2 2d2  60a2 þ   7 7 7 d ðd  2aÞ ðd þ 2aÞ d ð13Þ

where Acsc is the Hamaker constant and r is the Lennard-Jones diameter of atoms creating the NPs. Interaction between two NPs occurs across solvent (water), thus the Hamaker constant Acsc is a function of solvent properties, particularly its Hamaker constant Ass. By utilizing standard mixing rules for the Hamaker constants [17] we can easily find

Acsc ¼

pffiffiffiffiffiffiffi pffiffiffiffiffiffi2 Acc  Ass

ð14Þ

where Acc and Ass are, respectively, the Hamaker constants for interactions between colloid material and solvent material across vacuum. These values, determined experimentally, are available in literature for a number of materials [17]. They can also be calculated with a good accuracy from refractive indices of the considered materials using the Lifshitz theory of van der Waals interactions [17,18]:

Acc ¼

3 kT 4



2

ec  1 ec þ 1

2

þ

3hme ðn2c  1Þ pffiffiffi 16 2 ðn2c þ 1Þ3=2

ð15Þ

where ec is the dielectric constant of the material, nc is its refractive index, h is the Planck constant and me is absorption frequency normally taking the value 3  1015 s1. The Hamaker theory predicts also the following expression for the interaction energy of colloid NPs with solvent molecules

U CS ¼

2a3 r3 Acs 2 3

9ða2  d Þ

"

2

1

4

6

ð5a6 þ 45a4 d þ 63a2 d þ 15d Þr6 6

15ða  dÞ ða þ dÞ

6

# ð16Þ

pffiffiffiffiffiffiffiffiffiffiffiffi In this case Acs ¼ Acc Ass and solvent molecules are considered as simple point particles. In the case of core–shell architecture of the NPs the dispersion interactions across shell material are screened. Thus, the total interaction energy between two core–shell NPs depends on the

A colloidal suspension normally consists of macroions (charged NPs), counterions, salt ions, and solvent molecules. This multicomponent mixture can be modeled on an effective level within the DLVO theory [20–22]. It involves only the macroions, whose Coulomb repulsion is screened exponentially by the surrounding counterions and salt ions. The electrostatic interaction between two charged NPs is described by the screened Coulomb interaction given by the Yukawa term [22],

U Y ðrÞ ¼

ðZe0 Þ2 expð4jaÞ ejd 4pes e0 ð1 þ 2jaÞ2 d

ð17Þ

where Z is the number of elementary charges (e0) the macroion gains in the solution, e0 is the dielectric permittivity of vacuum, and es is the relative dielectric constant of the solvent. The inverse Debye screening length j involves contributions from the counterions and the salt ions. The salt contribution is measured by the ionic strength Is (in mol L1) [22]. Assuming charge neutrality between macroions and counterions, the inverse Debye length can be written as

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1000e20 j¼ ðZ q þ 2Is NA Þ e0 es kT

ð18Þ

where q is the macroion concentration (in mol L1) and NA is Avogadro’s number. Because the NPs are intended to work in biological environment we consider only the physiological ionic strength, Is = 0.145 mol L1. The amount of charge collected on NPs surfaces can be linked with a measurable quantity, that is zeta potential, f. Assuming that a NP creates a spherical capacitor which walls are separated by j1 then the charge collected on the walls becomes the following function of the voltage between walls, that is zeta potential:

Ze0 ¼ 4pe0 es að1 þ jaÞf

ð19Þ

Eq. (19) is obviously an approximation, however, more advanced approach [23] predicts only 3% larger value of the electroforetic charge at the considered conditions for f = 25 mV. 2.4. Proton transverse relaxation time The transverse relaxation time T2 of protons is affected by the presence of magnetic NPs in a sample due to magnetic field inhomogeneities. At high magnetic fields, the spin–spin relaxation is caused by local dephasing of each proton spin diffusing around magnetic NPs. In literature there are several papers devoted to determination of T2 time using the random walk method [6–10]. This method can be shortly summarized as follows. An ensemble of noninteracting motionless magnetic NPs is defined, they might be isolated from each other or form clusters of various shapes and sizes. Into that system some number of probe protons is inserted and they are allowed to diffuse according to the random walk scheme. That is, the position of each proton (formally belonging to water molecule) is updated at eachffi timestep Dt by adding a pffiffiffiffiffiffiffiffiffi randomly oriented vector of length 6DDt to the previous position. The diffusivity D corresponds to water molecules. The protons do not interact with each other and their interaction with NPs is only of hard sphere repulsion.

190

T. Panczyk et al. / Journal of Colloid and Interface Science 437 (2015) 187–196

During the diffusion, at high static external magnetic field aligned along the z-axis, the local field felt by protons spins induces a local Larmor precession which results in rotation of magnetic moments in the xy-plane, that is dephasing. This dephasing is proportional to the z-component of the local magnetic field felt by a given proton,



Bz ðr; hÞ ¼

~ j 3cos2 ðhÞ  1 l0 jm 4p r3

 ð20Þ

where r is the position relative to center of the NP and h is the angle between m and r. At high static fields all magnetizations are aligned with the field direction, thus h can be viewed as the inclination angle of r to the z-axis (field direction). The total Bz acting on a given proton, Btot, is the sum of contributions from all NPs. Thus, the dephasing for the i-th proton spin is given by

D/i ¼ cBtot ðiÞDt

ð21Þ

where Dt is the integration timestep and the cumulative phase at time t

/i ¼

t X D/i ; /i ð0Þ ¼ 0

ð22Þ

0

The normalized transverse magnetic moment of the i-th proton in the rotating frame is thus ~ li ¼ ðcos /i ; sin /i ; 0Þ and the normalized signal is then an average value of the transverse magnetic moment over all spins for each time t,

hlðtÞi ¼

1 Nspins

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi X 2 X 2 cos /i þ sin /i

ð23Þ

The mean transverse magnetic moment hli is then fitted by exp(t/T2) giving finally the simulated transverse relaxation time or the relaxation rate R2 ¼ T 1 2 . The above methodology was applied for analysis of T2 changes as a function of NPs sizes, distribution within the simulation box and their agglomeration into clusters of various sizes and shapes [7–10]. It was also verified by comparing its results with theoretical predictions available for some limiting cases. Thus, we use the same methodology for determination of T2 but we will focus on the factors which are not tractable using the random walk method. The random walk approach is very efficient in terms of timescale ranges possible to probe in simulations. By taking sufficiently large Dt the simulation can reach macroscopic time regimes. However, very large Dt means very large jumps of diffusing protons and it makes the assumption of independent motion questionable. Intermolecular interactions between NPs and water protons, as well as between NPs are not normally accounted for in the random walk method. Their meaning is difficult to predict without dedicated calculations. The approach utilized in this study is based on full molecular dynamics description of the system. It allows for coherent solution of motion equation according to Langevin dynamics and intermolecular forces acting in the system. The effects of magnetization reversal within the NPs and thermal fluctuations around the equilibrium directions can also be accounted for. 2.5. Simulation protocol The transverse relaxation times due to field inhomogeneity caused by the NPs are normally of the order of milliseconds. It means that we need to reach a comparable timescale in molecular dynamics calculations. Therefore, we need to impose several restrictions to the simulation box in order to make the problem tractable.

Thus, the calculations involve only several dozens of NPs (from 8 to 64) suspended in an implicit solvent. The calculations were based on the Langevin dynamics with the integration time step 0.01 ps. The damping factor for the probe protons was determined in such a way that the mean square displacement of the protons corresponded to water diffusivity, i.e. 2.3  105 cm2 s1. The damping factor for the NPs was calculated from the Einstein– Stokes equation. The implicit solvent parameters correspond to physiological fluid at standard density and ionic strength. The water protons, which measure the dephasing, are treated as the TIP3P water molecules moving in the implicit solvent and do not interact with each other. That simplification is allowed since the dephasing is only a function of the instantaneous position of probe protons. On the other hand, interactions between water molecules enter indirectly via diffusion coefficient. The number of probe water molecules was 128 or more depending on the number of the NPs. The NPs parameters varied depending on the factor being studied but normally we were close to the parameters’ ranges of iron oxide nanoparticles covered by ultrathin layer of modified chitosan. These NPs have been recently fabricated and carefully studied as promising contrast agents in MRI [3]. Thus, we assumed that individual NPs have diameters 120 Å and the thickness of chitosan layer is 10 Å. Their saturation magnetization was experimentally found to be 123 emu g1 Fe, which upon recalculation taking into account the density of magnetite gives Ms = 470 kA m1. The magnetic anisotropy constant was assumed to be 1.1  104 J m3, that is typical for bulk magnetite. The Hamaker constant for chitosan layer was calculated from Eq. (15) using as inputs the refractive index and dielectric constant found in the literature [24,25] and it takes the value Acc = 8  1020 J. The size of the simulation box (with periodic boundary conditions) varied depending on the number of the NPs in a given run and keeping the volume fraction of the NPs constant, that is, U = 0.01. Thus, for 8 NPs case the size of the box was l = 893 Å. In any case the cutoff for intermolecular interactions and for dephasing measurements was 0.5l. Such a long cutoff highly slows down the computations; however, as found in preliminary studies the cutoff of at least 5a was necessary for colloid–colloid interactions. Also, Bz decays slowly with the distance therefore so large cutoff was necessary. The calculations started from 10 ls equilibration runs and next 10–20 independent production runs were performed for any combination of system parameters. The final results were averaged over all independent runs thus the total simulation time for a given case was 0.1–0.2 ms. To sum up, the overall simulation scheme consists of two steps: (a) molecular dynamics simulation loop and (b) data analysis. Within the MD simulation loop (a) the system runs according to the assumed force field, i.e. the torques and forces acting on all species are calculated from Eqs. 3–5 and (13), (16) and (17), respectively. The rotational degrees of freedom of NPs and their magnetizations are integrated according to Eqs. (6), (7) and (12) whereas the coupling of the Néel and Brownian rotations enters the protocol via Eq. (11). The translational degrees of freedom are integrated according to the standard leapfrog algorithm. During the MD simulation loop (a) the following data is collected: temperature, radial and angular distribution functions, mean squared displacements and particularly local Bz for each proton Eq. (20), local dephasing Eq. (21) for each proton and its cumulative value Eq. (22). Every some number of integration timesteps the average transverse magnetic moment of the system is outputted according to Eq. (23). Next, beyond the MD loop (b) the data analysis is performed, that is the average transverse magnetic moment (Eq. (23)) is plotted as a function of time. Next, the fitting of the exponential function is performed and from that fit the time constant of the exponential decay is determined, that is the T2 value.

T. Panczyk et al. / Journal of Colloid and Interface Science 437 (2015) 187–196

191

3. Results and discussion 3.1. Validation of the model The applied simulation protocol uses a new computational methodology and therefore it needs validation before performing further analysis. Such a validation is possible by comparison of simulation results with exact solution when it is available. In the case of spherical noninteracting NPs the outer-sphere theory provides such a solution [6,9,10]. The outer-sphere theory is applicable when the magnetization of the NP is locked along the z-axis due to interaction with the external field, that is, the field strength must be high. Another necessary condition is the Redfield condition, that is

cBeq sD < 1

ð24Þ

where sD ¼ r 2M =D is water diffusion correlation time (D – diffusivity of water molecules) and

rffiffiffi 4 l0 Beq ¼ Ms 5 3

ð25Þ

is the equatorial field generated by a magnetic nanoparticle of radius rM. It is normally assumed that Beq for magnetite nanoparticle is about 0.16 T. The Redfield condition requires thus that the time necessary for diffusion of water molecule around the NP must be shorter than the time of proton magnetic moment rotation induced by the field generated by the NP. Thus, in principle, the outer-sphere theory is only applicable to very small nanoparticles. That limiting case is often called as motional averaging regime, MA, provided that no RF echo-pulses are applied or the time between pulses is much longer than sD. There are also other theories valid for larger NPs but we will not discuss them since we will limit our further analysis to the motional averaging regime. According to the outer-sphere theory the relaxation rate caused by the presence of the magnetic nanoparticle and within the MA regime reads,

1 16 ¼ UsD ðcBeq Þ2 T 2 45

ð26Þ

Thus, it possible to predict the relaxation time for a given set of parameters and compare it to the simulation results obtained in the same conditions. Therefore, we started from the analysis of the system which obeys the assumptions of the outer sphere model, that is, we consider a set of 8 NPs without protective shells so that the total NP volumes belong to magnetic material. For the purpose of validation tests we intentionally assumed that Hamaker constant for dispersion interactions is the same as if the NPs would be covered by chitosan layer (Acc = 8  1020 J) and they carry a charge corresponding to zeta potential found for chitosan covered NPs, that is 47 mV. This ensures that the ensemble of the NPs will not agglomerate in a cluster. The dispersion interactions between water molecules and NPs were intentionally reduced by 2 orders of magnitude (Acs was divided by 100) just to avoid adsorption of water on the NPs surfaces. These assumptions are obviously physically inconsistent but they allow mimicking the conditions of the outer sphere model in our simulation setup. Fig. 1 shows the results of validation tests, that is, the time dependence of the average value of the transverse magnetic moment hli and radial distribution functions for nanoparticle– nanoparticle (NP–NP) and nanoparticle–water (NP–w) proton distances. The hli vs. time was averaged over 20 independent runs from the initially equilibrated systems. As we can see, we have got a smooth exponential decay of hli vs. time. From fitting an

Fig. 1. Transverse magnetic moment as a function of time obtained from MD simulations for the following set of system parameters: rM = a = 60 Å, Ms = 470 kA m1, B = 7 T, Acc = 8  1020 J, Acs  0, f = 47 mV, U = 0.01. The T2 = 11.89 ls relaxation time is obtained from fitting exponential function to hli vs. time. The radial distribution functions, rdf, indicate that at the assumed conditions the density of NPs within the simulation box is uniform, that is, the NPs move independently without forming any clusters. The rdf for NP – water distances show that the probe water molecules do not adsorb on the NPs surfaces, that is, they are moving without retardation induced by dispersion interaction with the NPs. The inset shows a simulation snapshot; yellow spheres are the magnetic NPs and the red spots on their surfaces indicate the directions of their magnetizations while black spots (little visible as they are rotated by 180 deg) indicate the direction of the easy axes. Blue points show the positions of the probe water molecules. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

exponential function we obtained the T2 value for the assumed set of simulation parameters. That value, i.e. T2 = 11.89 ls is to be compared with T2 calculated from Eq. (26). By assuming Beq = 0.16 T and sD = 13.2 ns we find the T2 value predicted by the outer sphere model equal to 11.42 ls. The diffusion correlation time sD was calculated using the effective water diffusivity found in the simulations, that is, D = 2.73 ± 0.23  105 cm2 s1. That diffusivity differs slightly from the exact value for water due to application of the Einstein–Stokes equation for determination of the damping factor in the Langevin dynamics. The Redfield condition is satisfied at the studied conditions because the lhs of Eq. (24) is 0.56. Also, we do not observe any effects coming from intermolecular interactions because the densities of NPs as well as water molecules are uniformly distributed over the simulation box. Thus, the properties of the simulated system obey the assumptions of the outer sphere model. Thus, we can conclude that the simulation protocol correctly reproduces the theoretical prediction and can be used for further studies involving more complex situations. A slightly larger T2 than the outer-sphere prediction obtained from the simulations is a common effect found in other computational studies [10]. 3.2. Effects of intermolecular interactions In order to understand the influence of various intermolecular interactions parameters on the possible changes in relaxivity caused by the presence of magnetic NPs we studied several combinations of those parameters leading to distinct physico-chemical characteristics of the systems. To the best our knowledge some of them have never been analyzed in the literature, at least by using computer simulations. Table 1 collects the assumed simulation parameters for all case studies performed in this work. Let us first analyze the cases A–C which give us a notion about the influence of dispersion forces acting between NPs as well as between water protons and the NPs. Fig. 2 shows how the trans-

192

T. Panczyk et al. / Journal of Colloid and Interface Science 437 (2015) 187–196

Table 1 Definitions of the case studies performed in this work. If not pointed out strictly the following parameters values are the same for all analyzed systems: T = 300 K, is the transverse relaxation time predicted by the outer sphere model. Ms = 470 kA m1, f = 47 mV, U = 0.01, Ka = 1.1  104 J m3. T(OSM) 2

a

id.

nNP

nw

Acc (1020 J)

Acs (1020 J)

rM (Å)

a (Å)

B (T)

DNP (107 cm2 s1)

T2 (ls)

Comment

Test A B C D E F G H I J K L M N O P

8 8 8 8 8 8 8 8 8 8 8 4 8 18 27 48 64

128 128 128 128 128 128 128 128 128 128 128 128 128 288 432 768 1024

8 8 8 21 21 8 8 8 21 8 8 8 8 8 8 8 8

0 11 11 18 0 11 11 11 18 0 11 11 11 11 11 11 11

60 60 50 60 60 50 60 50 60 60 50 50 50 50 50 50 50

60 60 60 60 60 50 60 60 60 60 60 60 60 60 60 60 60

7 7 7 7 7 7 0.5 0.5 0.5 0.5 0.3 7 7 7 7 7 7

6.16 ± 1.46 4.3a – 7.52 ± 0.88 1.50 ± 0.14 – 7.19 ± 0.23 – – – – – 2.86 ± 1.10 1.44 ± 0.59 0.59 ± 0.20 0.37 ± 0.14 0.21 ± 0.17 0.15 ± 0.04

11.89 ± 1.34 11.06 ± 1.34 34.77 ± 3.51 6.49 ± 1.10 5.26 ± 0.77 16.20 ± 1.17 11.89 ± 1.60 33.84 ± 4.38 6.52 ± 0.88 11.63 ± 1.32 39.28 ± 4.55 18.11 ± 2.63 12.49 ± 1.66 11.26 ± 0.92 10.58 ± 1.25 11.82 ± 0.93 12.06 ± 0.74

T(OSM) = 11.42 ls 2 Free NPs Free NPs Agglomeration Agglomeration Free NPs Free NPs Free NPs Agglomeration Free NPs Free NPs Rigid cluster Rigid cluster Rigid cluster Rigid cluster Rigid cluster Rigid cluster

Calculated from the Einstein–Stokes equation, D ¼ 6pkTga, g dynamic viscosity of water.

Fig. 2. Transverse magnetic moments and radial distribution functions (rdf) determined for the systems A–C and using the parameters collected in Table 1. The code ‘A/N’ means the rdf for the NP–NP distances for the system A whereas ‘ A/w’ means the rdf for the NP–water molecule distances for the system A.

verse magnetic moment of protons decays with time in these cases, it also shows structural characteristics of the systems, i.e. the radial distribution functions. The case A corresponds to quite artificial situation because the NPs reveal the Hamaker constant representative of chitosan (or

other material with similar Acc) and saturation magnetization of iron oxide. Thus, the whole volume of the NP belongs to magnetic material, rM = a, but the dispersion interactions are screened as it would be covered by chitosan. However, it is interesting to compare the system A with the test system corresponding to the outer-sphere model. The difference between these two systems concerns only the Hamaker constant for NP–water interactions (nonzero Acs determined from the standard mixing rule pffiffiffiffiffiffiffiffiffiffiffiffi Acs ¼ Acc Ass ). Thus, in the case A the probe protons carried by water molecules are able to adsorb on the NPs surfaces. Indeed, there appears a density excess at the vicinity of the NPs surfaces in the case A, whereas in the test case the rdf was completely flat. The adsorption is, however, fully reversible as observed in the simulations snapshots. The density excess is due to retardation of water motion being the result of the attraction between NPs and water molecules. Because in actual measurements the adsorption of water is a common phenomenon its effect in modification of the relaxation rates needs to be discussed. As we can see the adsorption of water reduces T2 (enhances relaxation rate) and this is fully understandable because the probe protons in such case stay at the vicinity of the NPs for longer time and dephasing occurs at higher local fields. However, the difference between T2’s is surprisingly small. In fact, it is within the estimation errors so the adsorption of water has negligible influence on the relaxation times. For stronger interactions, that is, for cases C and D the difference is slightly higher but still not big and probably without a meaningful effect in experiments. It should be noted that cases C and D correspond to huge dispersion interactions between NPs and between water molecules and NPs. It leads to agglomeration of the NPs into a single cluster which does not decompose within the studied simulation times. The Hamaker constant for C and D corresponds to iron oxide [17] thus these cases are representative for bare magnetite nanoparticles without any protective layer. As we can see in the rdf in Fig. 2 the distances between NPs take discrete values with a very sharp maximum at 120 Å (closest contact distance) and much smaller at 200 Å. The latter value is very close to the second nearpffiffiffiffiffiffiffiffi est neighbor distance in bcc cubic crystal, that is 4 2=3a ¼ 196A. Thus, we observe agglomeration of the NPs and because any other parameter (except of Acc) is the same as in the case A the observed strong reduction of T2 must be attributed to the presence of cluster instead of freely moving isolated NPs. The reduction of T2 as a result of agglomeration of NPs into a cluster, while keeping the density constant, is well known in the literature [9,7,10].

T. Panczyk et al. / Journal of Colloid and Interface Science 437 (2015) 187–196

193

The case B corresponds to the situation when 100 Å in diameter magnetic core is covered by 10 Å thick chitosan layer. Thus, this is the case equivalent to the NPs studied in Ref. [3] In this case we observe much longer T2 than in cases A, C and D. Obviously this is the effect of reduction of the magnetic core size and steric barriers preventing close contacts of water protons with the magnetic cores surfaces. By comparing the case B with E we can notice a significant reduction of the relaxation rate due to the presence of the nonmagnetic protective layer. Again, the case E is physically artificial because the NPs interact as they would be covered by chitosan but the nonmagnetic steric barriers are removed. Thus, we can conclude that chitosan (or other nonmagnetic) layer must be as thin as possible because it strongly enhances T2. This is normally not proffered in a design and fabrication of the MRI contrast agents. 3.3. Effect of the field strength The field strengths applied in MRI cover quite wide ranges of values. Generally, the MRI is divided into high-field and low-field MRI. Both solutions have advantages and disadvantages, however, a detailed discussion on this matter is beyond the scope of this work [26]. The described simulation protocol allows for direct analysis of the magnetic field strength in modification of the transverse relaxation rates. That factor is not normally accounted for as it is believed that typical field strengths in MRI are enough to keep magnetizations fixed along the field direction. So, in order to verify such an assumption we determined T2 relaxation times at various field strengths but keeping all other parameters the same. Most calculations performed in this study assume high field limit >1.5 T, that is, 7 T. The cases F–J correspond to lower fields, that is, 0.5 T and 0.3 T (low-field MRI). The weakest field i.e. 0.3 T, approaches the limit of applicability of the assumed simulation protocol. However, the term cBsD is still slightly greater than unity at 0.3 T (1.06). Thus, according to the condition provided in [10] the protocol should still be reliable at 0.3 T. The T2 values for cases F–J (Table 1) should be compared to their high-field counterparts A–D. That comparison leads to the conclusion that the relaxation rates are not significantly altered when going from high to low (0.5 T) field MRI. The differences between T2 in both B regimes are within the estimation errors. However, for still lower fields, that is 0.3 T, the T2 strongly increases. This is a result of thermal fluctuations of magnetizations orientations of individual NPs. Fig. 3 shows how the interaction energy with the ! ~  B i changes with the field strength. For large fields that field hm energy is very high and this leads to very sharp distribution of ~ j Þ with the maximum for parallel alignment. ~ i; m the angles u ¼ \ðm For 0.3 T the energy drops to 44 kJ mol1 and the angle distribution function becomes more diffuse. This means more frequent small displacements of magnetization orientations from the field direction and thus the local Bz felt by water protons is significantly reduced. The mean magnetic anisotropy energy hEai shown in Fig. 3 is very small in any case. It is below kT for the case of the studied NPs. Thus, the meaning of the magnetic anisotropy is actually negligible for small NPs. An interesting feature is, however, the decrease of hEai for 0.3 T. This is due to the increasing meaning of rotational diffusion at low fields. The Brownian rotation is additionally driven by stochastic forces coming from collisions with solvent particles. In turn, it affects the Néel rotation due to the coupling in Eq. (11). 3.4. Effect of cluster size As found in Section 3.2, bare magnetite NPs undergo agglomeration into clusters at considered conditions. At the same time, we

Fig. 3. (Top): angle distribution functions for magnetizations orientations and the corresponding T2 dependence on the field strength. (Bottom): mean magnetic ! ~  B i. anisotropy energy hEai and mean energy of magnetization–field interaction hm

observed strong reduction of T2 and that effect must be attributed to the presence of cluster instead of freely moving isolated NPs (case C vs. A). As this is important observation we performed more detailed analysis of that phenomenon. Fig. 4 shows a series of clusters generated for that purpose. The clusters were generated by assuming that initially the NPs are uncovered by protective layers and reveal the Hamaker constant representative to magnetite, that is, 21  1020 J and their whole volumes belong to magnetite, so a = rM = 60 Å. Next, the NPs were allowed to relax until they finally agglomerated due to strong dis-

Fig. 4. Snapshots of rigid clusters analyzed in case studies K–P. The numbers at the bottom of each snapshot show how many individual NPs contribute to a given cluster whereas the numbers in square brackets show their diameters determined as 2RG, where RG is the gyration radius.

194

T. Panczyk et al. / Journal of Colloid and Interface Science 437 (2015) 187–196

persion forces. To keep the density constant the size of the simulation box for various numbers of the individual NPs was suitably increased. Also the number of probe protons was accordingly increased to keep their density the same in all case studies. Actually it was not necessary because the mean transverse magnetic moment was independent of the number of water protons, however, with larger number of protons the resulting hli vs. t was smoother though strongly enhanced computation time. Agglomeration process took about 1 ls for the largest system and afterwards the production runs were started. In the production runs the sizes of the magnetic cores were reduced to rM = 50 Å while the total radii a were kept equal to 60 Å. Thus, it means that the NPs received 10 Å thick chitosan layers with Acc = 8  1020 J. At those conditions the clusters would decompose because systems parameters correspond to the case A. Thus, in order to keep the clusters intact were applied the rigid body motion protocol to the clusters. It means that all distances, angles and interactions between NPs forming the rigid body were frozen and their equations of motions were reduced to single set corresponding to motion of a single body with its own inertia. In that way, we obtained freely moving clusters of chitosan covered magnetite NPs with frozen distances and interactions between NPs creating a given cluster. Fig. 5 shows how T2 changes with the cluster size or with the number of NPs creating the cluster. As we can see there exists a distinct minimum of T2 for cluster composed of 27 NPs, that is having the diameter about 33 nm. This is in good agreement with the results published by Vuong et al. [10] concerning the R2 dependence on the NPs radii. They found maximum of R2 (minimum of T2) for nanoparticles having radii about 20 nm and those NPs sizes correspond to an intermediate region between the motional averaging regime and partial refocusing one [10]. As seen in Fig. 5 we

observe similar effect for clusters of NPs which are additionally covered by nonmagnetic protective layers. Thus, we may conclude that there exists an optimal cluster size for which T2 is minimized for a given set of NP parameters and this optimal size (diameter) is in the range 30–40 nm for the considered magnetite particles. The decrease of T2 when going from single particles to small clusters is most pronounced as we observe more than 3 times larger T2 for single NPs compared to the minimal value for the cluster N27. Further increasing of the cluster size enhances T2 but this not very strong effect, at least at the considered clusters sizes. On the other hand, large clusters reveal small values of diffusivities and this obviously affect their mobility and perhaps biodistribution when thinking about their use as contrast agents. The bottom part of Fig. 5 shows the log–log plot of T1 2 vs. NNP. As can be seen there is a power law correlation between the relaxivity and the number of particles creating the clusters for small clusters consisting of about 27 NPs. That power law scaling has been found by Brown et al. [9] in their Monte Carlo studies focused on the aggregation of NPs and their influence on the relaxation rate. They also found that in case of clusters formed by diffusion limited aggregation (actually the same mechanism like in our case) the mean value of the power law exponent within the motional 0.44 averaging regime is 0.44, that is T1 . The dashed line in 2  NNP Fig. 5 shows that correlation in the case of our data. Clearly, the power law scaling with the exponent 0.44 holds in our case. The deviations observed for big clusters can be related to the dependence of T2 on cluster shape, as found by Brown et al. [9]. That effect is more pronounced in the case of big clusters, thus, in order to get a perfect power law scaling in the whole region we would have to collect a large set of T2 for the same NNP but for various clusters shapes and take its mean value. 3.5. Analysis of experimental data

Fig. 5. Transverse relaxation time T2 and diffusivity DNP for clusters (top) composed of NNP individual nanoparticles. Power law scaling of T1 with the number of 2 nanoparticles in the cluster (bottom).

Chitosan covered iron oxide NPs have been recently fabricated and characterized by means of several experimental techniques [3]. Particularly their r2 relaxivities were determined in order to assess their potential applicability as MRI contrast agents. In this study we tried to keep the parameters of the NPs close to the SPION particles analyzed in Ref. [3] and thus we are able to draw some conclusions concerning their properties from the microscopic point of view. The SPION particles were obtained by co-precipitation of ferrous and ferric salt solutions in the presence of cationic derivative of chitosan. In that way positively charged particles revealing +47 mV zeta potential were fabricated. In the second step those particles were treated by anionic derivative of chitosan using layer-by-layer (LbL) deposition technique. In this way negatively charged SPION were obtained revealing zeta potential 41 mV. In both cases the thicknesses of chitosan layers were very small, probably not larger than 10–30 Å, which is typical for LbL films formed as a result of electrostatic self-assembly of polyelectrolytes. The average diameters of particles in both cases were similar, i.e. 10–12 nm. Considering the above information we can notice that the system B in Table 1 reveals the same or similar parameters as the SPION nanoparticles. As indicated in Table 1 the NPs in case B move freely and do not undergo agglomeration. This is because the thin chitosan layers effectively screen the dispersion forces between magnetic cores and there is additional electrostatic repulsion coming from the equal sign charges. Because the TEM images of SPION show that they exist in the form of small aggregates we can conclude that creation of aggregates must occur during co-precipitation. Alternatively, the bare iron oxide NPs created during co-precipitation may form large aggregates which, upon addition of solution of chitosan

T. Panczyk et al. / Journal of Colloid and Interface Science 437 (2015) 187–196

derivative, dissolve due to adsorption of chitosan forming small clusters covered by chitosan layer. Finally, the resulting clusters are highly stable in aqueous solution due to the screening of dispersion interactions. At the same time the interactions between the magnetic cores within the clusters are very strong and they keep the internal structure of clusters intact. We can thus notice that the above mechanism of aggregation is equivalent to the case studies K–P. For both cases of SPION nanoparticles the experimentally determined r2 relaxivities were very high; 227 mM1 s1 for positively charged particles and 369 mM1 s1 for the negatively charged ones. Recalculation of T2 times obtained from simulations into r2 values available form experiments is not straightforward because it needs detailed knowledge about iron content in the NPs. However, just to have a notion about r2 ranges corresponding to the determined T2 values we assumed the stoichiometry of Fe3O4 for magnetic cores and the densities 5.3 g cm3 and 0.3 g cm3 for magnetite and chitosan layer, respectively. Table 2 shows r2 values calculated under the above assumptions for selected cases listed in Table 1. As we can see, the r2 values compare well with the experimentally determined ones though it seems that all of them are slightly lower than the reported experimental values, particularly for the negatively charged SPION. This is probably due to slightly smaller diameters of magnetic cores assumed in simulations or to the overestimation of iron content while recalculating T2 into r2. Nevertheless, it is clear that results obtained from simulations, which correspond to large volume fractions U = 0.01 of NPs, are reliable and can be extrapolated to realistic U values of the order 106. However, more important than a perfect match between experimental and theoretical values is the trend in r2 values as a function of the NPs diameters. As we can see the maximum in r2 corresponds to the clusters of sizes about 300–350 Å. The experimentally studied SPIONs reveal hydrodynamic diameters about 1000 Å [3], thus even if we account for the presence of hydration layer we still observe that they are definitely larger than the optimal cluster size found in simulations. So, we can conclude that the performance of SPION nanoparticles as MRI contrast agents could still be improved by modification of their synthesis method leading to production of smaller clusters with diameters about 30–40 nm. The difference between r2 for positively and negatively charged SPIONs cannot be directly explained using the simulation results. As mentioned already, the sign of the NP charge is not important, also some small difference in absolute values of zeta potentials cannot explain such large differences in r2 values. The most reasonable explanation is that the clusters differ in magnetic cores densities. The negatively charged SPIONs probably reveal denser packing of iron oxide nanoparticles than the positively charged ones. As a result, for those denser clusters the magnetic flux density Bz is higher and produces faster decay of transverse magnetic moment of protons. This explanation seems to be supported by

Table 2 Relaxivities r2 determined for the selected case studies shown in Table 1 by assuming the stoichiometry of Fe3O4 for magnetic cores and the densities 5.3 g cm3 and 0.3 g cm3 for magnetite and chitosan layer, respectively. The diameters were determined from the gyration radii of clusters. id

Diameter (Å)

r2 (mM1 s1)

B C K L M N O P

120 205 169 205 303 335 441 503

71 382 137 198 218 239 204 201

195

TEM images of both clusters reported in [3]. The positively charged SPION cluster is more loosely when compared to its negatively charged counterpart. Also, analysis of TEM images leads to the conclusion that an average number of individual NPs included in a single cluster is greater in the case of negatively charged SPION (probably a few hundreds) than in the case of positively charged ones (probably less than one hundred). 4. Summary In this work we proposed a method of determination of transverse relaxation times of protons based on molecular dynamics simulations. The method allows for studying such factors as dispersion interactions between NPs and water protons, magnetic field strength and thermal fluctuations of magnetic moments being due to finite magnetic anisotropy constants as well as core–shell architecture of the NPs. Generally, the model represents a full molecular dynamics force field representative of the suspension of magnetic nanoparticles in aqueous solution. We found that dispersion interactions between the NPs and water protons affect T2 relaxation times but that effect is rather small and without practical importance. On the other hand, the dispersion interactions between NPs are highly important and, depending whether the NPs are covered by a protective layer (with a moderate Hamaker constant) or bare, they might move freely in the solution or form agglomerates. The agglomeration, in turn, strongly affects T2 relaxation times. We found that there exists an optimal cluster size (300–350 Å) for which T2 is minimized. Analysis of the external magnetic field strength led to the conclusion that in typical MRI applications the thermal fluctuations of magnetic moments are not significant, in turn, they might be important in low field MRI because thermal fluctuations strongly enhance T2 relaxation times at low magnetic fields. Acknowledgments This research has received funding from the Marian Smoluchowski Krakow Research Consortium – a Leading National Research Centre KNOW supported by the Ministry of Science and Higher Education. L. K. and T. P. thank the National Science Centre (NCN) for financial support within the grant UMO-2012/07/E/ST4/ 00763. This work was also supported by the European Union from the resources of the European Regional Development Fund under the Innovative Economy Programme (grant coordinated by JCETUJ, No POIG.01.01.02-00-069/09). References [1] C.F.G.C. Geraldes, S. Laurent, Contrast Media Mol. Imaging 4 (2009) 1–23. [2] S. Laurent, D. Forge, M. Port, A. Roch, C. Robic, L. Vander Elst, et al., Chem. Rev. 108 (2008) 2064–2110. [3] A. Szpak, G. Kania, T. Skórka, W. Tokarz, S. Zapotoczny, M. Nowakowska, J. Nanoparticle Res. 15 (2012). [4] R.A. Brooks, F. Moiny, P. Gillis, Magn. Reson. Med. 45 (2001) 1014–1020. [5] R.A. Brooks, Magn. Reson. Med. 47 (2002) 388–391. [6] P. Gillis, F. Moiny, R.A. Brooks, Magn. Reson. Med. 47 (2002) 257–263. [7] Y. Matsumoto, A. Jasanoff, Magn. Reson. Imaging 26 (2008) 994–998. [8] N.P. Blockley, L. Jiang, A.G. Gardener, C.N. Ludman, S.T. Francis, P.A. Gowland, Magn. Reson. Med. 60 (2008) 1313–1320. [9] K.A. Brown, C.C. Vassiliou, D. Issadore, J. Berezovsky, M.J. Cima, R.M. Westervelt, J. Magn. Magn. Mater. 322 (2010) 3122–3126. [10] Q.L. Vuong, P. Gillis, Y. Gossuin, J. Magn. Reson. 212 (2011) 139–148. [11] S. Plimpton, J. Comput. Phys. 117 (1995) 1–19. [12] M.F. Hansen, S. Mørup, J. Magn. Magn. Mater. 184 (1998) L262–274. [13] J. García-Palacios, F. Lázaro, Phys. Rev. B 58 (1998) 14937–14958. [14] W. Scholz, T. Schrefl, J. Fidler, J. Magn. Magn. Mater. 233 (2001) 296–304. [15] W. Wernsdorfer, E. Orozco, K. Hasselbach, A. Benoit, B. Barbara, N. Demoncy, et al., Phys. Rev. Lett. 78 (1997) 1791–1794. [16] R. Everaers, M. Ejtehadi, Phys. Rev. E 67 (2003) 041710. [17] J.N. Israelachvili, Intermolecular and Surface Forces, Academic Press, Burlington, MA, 2011.

196

T. Panczyk et al. / Journal of Colloid and Interface Science 437 (2015) 187–196

[18] T. Panczyk, W. Rudzinski, A. Jagusiak, Colloids Surf. Physicochem. Eng. Asp. 409 (2012) 149–158. [19] T. Panczyk, P.J. Camp, G. Pastorin, T.P. Warzocha, J. Phys. Chem. C 115 (2011) 19074–19083. [20] J.-P. Hansen, Theory of Simple Liquids, third ed., Elsevier/Academic Press, Amsterdam; Boston, 2007. [21] B. Fazelabdolabadi, J.Y. Walz, P.R. Van Tassel, J. Phys. Chem. B 113 (2009) 13860–13865.

[22] S. Grandner, S.H.L. Klapp, J. Chem. Phys. 129 (2008) 244703. [23] M. Tirado-Miranda, C. Haro-Pérez, M. Quesada-Pérez, J. Callejas-Fernández, R. Hidalgo-Álvarez, J. Colloid Interface Sci. 263 (2003) 74–79. [24] D.E. Azofeifa, H.J. Arguedas, W.E. Vargas, Opt. Mater. 35 (2012) 175–183. [25] Y. Li, W. Gao, K. Zhao, G. Yang, Z. Zhu, R. Cui, J. Phys. Chem. B 117 (2013) 3337– 3344. [26] M. Espy, A. Matlashov, P. Volegov, J. Magn. Reson. 229 (2013) 127–141.