View Online / Journal Homepage / Table of Contents for this issue

PAPER www.rsc.org/pccp | Physical Chemistry Chemical Physics

Heat transfer in protein–water interfaces

Anders Lervik,ab Fernando Bresme,*ac Signe Kjelstrup,bc Dick Bedeauxbc and
J. Miguel Rubicd
Received 8th September 2009, Accepted 3rd December 2009
First published as an Advance Article on the web 11th January 2010
DOI: 10.1039/b918607g

We investigate using transient non-equilibrum molecular dynamics simulation the temperature

relaxation process of three structurally different proteins in water, namely; myoglobin, green
fluorescence protein (GFP) and two conformations of the Ca2+-ATPase protein. By modeling the
Published on 11 January 2010 on http://pubs.rsc.org | doi:10.1039/B918607G

temperature relaxation process using the solution of the heat diffusion equation we compute the
thermal conductivity and thermal diffusivity of the proteins, as well as the thermal conductance of
the protein–water interface. Our results indicate that the temperature relaxation of the protein can
Downloaded by University of Memphis on 08 October 2012

be described using a macroscopic approach. The protein–water interface has a thermal

conductance of the order of 100–270 MW K1 m2, characteristic of water–hydrophilic interfaces.
The thermal conductivity of the proteins is of the order of 0.1–0.2 W K1 m1 as compared
with E0.6 W K1 m1 for water, suggesting that these proteins can develop temperature
gradients within the biomolecular structures that are larger than those of aqueous solutions. We
find that the thermal diffusivity of the transmembrane protein, Ca2+-ATPase is about three times
larger than that of myoglobin or GFP. Our simulation shows that the Kapitza length of these
structurally different proteins is of the order of 1 nm, showing that the protein–water interface
should play a major role in defining the thermal relaxation of biomolecules.

Introduction transfer between these sites may involve the concerted

motion of hundreds to thousands of atoms. Studies on
Biological molecules have developed mechanisms and struc- different proteins suggest that there might be a set of residues
tures that enable the efficient and rapid dissipation of excess that configure the energy pathway. This energy pathway
energy through the biomolecule structure itself and through would be an inherent property of the protein’s tertiary
the biomolecule–solution interface. Finding the microscopic
structure, and could provide the microscopic basis for signal
mechanism controlling energy transport in biomolecules
transduction.5
represents a considerable challenge, which could provide a
Chemical and photochemical reactions occurring at specific
basis to understand complex biological processes. It is accepted
spots in biomolecules can result in large increases in temperature
that the timescale associated with the energy relaxation
in very small volumes. Transient photon experiments, where
throughout the molecule can influence the kinetics of biomolecule
photon absorption is converted to vibrational energy, show
reactions.1 Similarly the pathways followed by the energy
that the temperature rise resulting from this process can be
during the relaxation process can provide important clues to
very significant, between 500 and 1000 K.6,7 Recent work on
understand allostery.2 In order to rationalize these processes,
the Ca2+-ATPase embedded in the sarcoplasmic reticulum has
some questions must be addressed, e.g., how perturbations
suggested that this enzyme can—under working conditions—
occurring at a specific spot in the biomolecule, namely, binding
release significant amounts of heat. Using micro-thermometers
of small molecules at receptor sites, may impart a conforma-
tional change at a distant spot, lying several nanometres it was found that thermal gradients of the order of 105 K m1
away.3 This question is very relevant to understanding can develop between regions separated by tens of micrometres.8
how biomolecular motors, such as Ca2+-ATPase, work. These temperature gradients could lead to interesting coupling
Ca2+-ATPase requires the hydrolysis of ATP to enable effects, such as water polarization, which has been predicted
Ca2+ transport. In this instance, the spot where the hydrolysis very recently.9
takes place and the binding site for the Ca2+ ions, are The quantification of thermal transport in biomolecules in
separated by about 4 nm.4 It is expected that the energy terms of a few bio-material properties, e.g. the thermal diffusivity,
has been discussed recently.10 Experimental and theoretical
a
studies have provided important insights on the vibrational
Department of Chemistry, Imperial College London, London,
energy relaxation and energy pathways in proteins. Using
UK SW7 2AZ. E-mail: f.bresme@imperial.ac.uk
b
Department of Chemistry, Norwegian University of Science and picosecond anti-Stokes resonance Raman spectroscopy it has
Technology, NO-7491, Trondheim, Norway been possible to monitor the cooling dynamics of myoglobin.11 It
c
Center for Advanced Study at the Norwegian Academy of Sciences is expected that in addition to the energy relaxation occurring
and Letters, NO-0271, Oslo Norway
d
Departamento de Fı´sica Fonamental, Facultad de Fı´sica, Universidad inside the protein, vibrational energy can also flow through the
de Barcelona, E-08028, Spain. E-mail: fbresme@imperial.ac.uk protein–solvent interface. Experiments of peptide chains in an

1610 | Phys. Chem. Chem. Phys., 2010, 12, 1610–1617 This journal is 
c the Owner Societies 2010
 View Online

organic solvent indicate that a significant amount of dissipation ATPase was the same as that of the vesicle interior. This
can occur through the molecule–solvent interface.7 The time assumption remains unsupported. Hence, in order to further
scale associated with such processes has also been estimated, establish the validity of these assumptions it is necessary to
being of the order of tens of picoseconds in myoglobin.12 This investigate the thermal properties of proteins, and study their
time scale is not far from the one associated to the relaxation dissipation of energy.
through the biomolecular structure, indicating that a vibrational The paper is structured as follows. Firstly we discuss the
coupling between biomolecular and solvent degrees of freedom methodology employed in this work; the transient non-
may take place.13 equilibrium molecular dynamics simulations and the equations
Computer simulation techniques are being used to investigate employed to analyze the temperature relaxation in the proteins.
the transport properties as well as vibrational energy pathways The discussion of our results including the heat capacity of the
of proteins.6,7,10,12,14–28 These studies have reported vibra- proteins as well as the thermal conductivity and thermal
tional energy relaxation times between one to several tens of conductance follows. We finish the paper with the main
picoseconds. The first simulation study of vibrational energy conclusions obtained in this work.
Published on 11 January 2010 on http://pubs.rsc.org | doi:10.1039/B918607G

relaxation of heme proteins was reported in ref. 6. These

simulations relied on the injection of vibrational energy by
Methodology
an amount equivalent to the energy associated to photon
absorption. Subsequent methodologies have exploited the Transient non-equilibrium molecular dynamics
Downloaded by University of Memphis on 08 October 2012

scaling of the energy diffusion coefficient with the vibrational

The interfacial conductance can be computed using stationary
mode frequency of a protein,10 or have modelled the heat
non-equilibrium molecular dynamics (NEMD) simulations,
transfer process as a boundary thermal transport problem.14–16
whereby a temperature gradient is imposed in a direction
These simulation methods offer the possibility of estimating
perpendicular to the interface plane.34–38 This methodology
transport coefficients, particularly the thermal diffusivity and
has been used to compute the interfacial conductance, or the
the thermal conductivity. In this way it has been shown that
inverse of this quantity—the interfacial resistivity—of alkane–
the thermal conductivity of several proteins, myoglobin,
water and alkane–vapor interfaces.37,38 In addition to stationary
ribonuclease T1, and green fluorescent protein GFP, is
methods, non stationary NEMD or transient NEMD has been
smaller than that of water, about 200–300 mW K1 m1.10,16,23
employed to investigate biomolecules.14,16 This method
This value is similar to the thermal conductivity of
explicitly models the protein–water interface by immersing
hydrocarbon liquids29 and hydrocarbon nanodroplets.30
the protein in a solvent. The protein and the solvent are
In addition to quantifying transport coefficients, computer
set initially at different temperatures, and the temperature
simulations provide a unique approach to interpret the
relaxation is monitored over time. The temperature relaxation
vibrational relaxation process in terms of intermolecular inter-
curve is fitted to a solution of the heat diffusion equation,
actions. As a matter of fact, it has been shown that the
which can be used to estimate the interface thermal conduc-
electrostatic interaction of the isopropionate side chains of
tance as well as the thermal conductivity and diffusivity. Many
myoglobin with the surrounding solvent, provides a mechanism
computational studies assume a first order process in which the
for kinetic energy dissipation.19,20
temperature relaxation is modeled as a simple exponential
Previous experimental and simulation studies indicate
decay. We have considered recently an analysis that goes
that the solvent surrounding a protein can enhance the
beyond a pure first order process, and that incorporates as a
dissipation of vibrational energy associated with high
variable the thermal conductivity and the temperature
temperature regions inside the protein. It is expected that the
discontinuity arising at the protein–water interface, due to
thermal conductance of the protein–water interface may
the acoustic mismatch of the protein and the solvent.30 This is
play an important role in determining the efficiency of this
the approach we use in this work. For completeness it is
relaxation process. We have recently shown30 that the inter-
discussed below.
face curvature strongly affects the magnitude of the interfacial
The temperature relaxation can be modeled using the heat
thermal conductance. This effect is particularly important at
diffusion equation (in spherical coordinates),
the nanometre scale characteristic of nanoparticles, and is  
hence expected to be relevant also in biological nanoparticles 1 @ 2 @Tðr; tÞ @Tðr; tÞ
kr ¼ rcp ð1Þ
(proteins). r2 @r @r @t
In this article we have used transient non-equilibrium
molecular dynamics to compute the thermal conductance of along with the boundary condition equation,
the protein–water interface for three structurally different 
@T 
k ¼ GðTðr ¼ R; tÞ  Tf Þ; ð2Þ
proteins, namely: myoglobin, the green fluorescence protein @r r¼R
and two conformations of the Ca2+-ATPase. The latter is an
example of an important family of molecular motors. Using where k and G are the thermal conductivity and thermal
the non-equilibrium thermodynamics theory, non-linear conductance, respectively, cp is the isobaric heat capacity, R
flux-force relations for conversion of chemical to osmotic is the radius of gyration of the protein, T(r, t) is the local
and thermal energy have been derived recently.31–33 It has temperature at position r and time t inside the protein, and Tf
been possible to obtain an expression for the measurable heat is the temperature of the surrounding fluid, which it is assumed
flux out of the ATPase and the vesicle it was embedded in. In to be constant. We note that the boundary condition given in
order to do this it was assumed that the temperature of the eqn (2) is equivalent to a convective boundary condition.

This journal is 
c the Owner Societies 2010 Phys. Chem. Chem. Phys., 2010, 12, 1610–1617 | 1611
 View Online

Using the initial condition, T(r, t = 0) = Ti, the solution of estimates for A and B. l1 is obtained from A and t from B. The
eqn (2) is39 thermal conductivity follows from
Tðr;tÞ  Tf
Ti  Tf rcp R2 3Cp
k¼ ¼ ð9Þ
 2   ð3Þ t 4pRt
X
1
sinðln Þ  ln cosðln Þ l t R ln r
¼ 4 exp  n sin ;
n¼1
2ln  sinð2ln Þ t ln r R where R is the radius of gyration of the protein, r is the protein
2 density, which we define as r = m/V, where m is the protein
where t = R /DT, DT is the thermal diffusivity, DT = k/(rcp)
mass and V = 4pR3/3, Cp = mcp (in J K1). Eqn (9) follows
The coefficients, ln, fulfil the following relation,
from DT = R2/t = k/(rcp). Once k is known, the thermal
GR conductance (G) can be estimated using eqn (7). In addition to
1  ln cotðln Þ ¼ ¼ Bi; ð4Þ
k the analysis of the temperature relaxation, the computation
of the thermal conductivity requires the heat capacity, Cp, and
Published on 11 January 2010 on http://pubs.rsc.org | doi:10.1039/B918607G

where Bi is the Biot number. The Biot number measures the

ratio of the thermal resistance of the protein to the thermal the radius of gyration of the proteins (cf. eqn (9)). The
resistance of the interface. A first order process normally computation of the heat capacity is discussed below.
corresponds to Bi { 1. Considering a thermal conductance Using the approach outlined above we have estimated the
Downloaded by University of Memphis on 08 October 2012

of G E 102 MW K1 m2, typical of hydrocarbon–water thermal conductance of small alkane nanodroplets in water.30
interfaces,30,37 a protein radius of gyration of E2 nm and a Our results for this quantity when extrapolated to a droplet of
protein thermal conductivity of the order of 0.2 W K1 m1,23 infinite radius were consistent with the estimates obtained for
we get Bi E 1. As we will see below, the thermal conductance planar alkane–water interfaces using the stationary NEMD
of the protein–water interface is similar to that of the approach.37 This lends support to the consistency of our
alkane–water interface. Hence, our estimate of Bi E 1 indicates method.
that the thermal resistance of the protein and the protein– We note that the equations discussed above (3–6) are used
water interface are of the same order, and both the thermal to model heat relaxation in macroscopic objects too. We find
conductivity and the thermal conductance should be included that the relaxation in small objects such as proteins can be
in the diffusion equation to model the temperature relaxation described in the same way as in larger macroscopic objects.
process. One requirement for this to be true is that the velocity
The solution of eqn (3) can be simplified by taking the components of the atoms must follow a normal distribution,
average temperature over the whole volume, V, of the protein. so that we can define a proper microscopic temperature.
This results in Deviations from normality if present are expected to appear
more clearly in very small objects. Hence, to test this point we
hTðr; tÞiV  Tf TðtÞ  Ti have chosen the smallest protein investigated in this work,
¼
Ti  Tf Ti  Tf myoglobin. Fig. 1 shows the normalized probability velocity
X1  2
1 ðsinðln Þ  ln cosðln ÞÞ2 ln t distribution as well as the quantile–quantile (Q–Q) plot
¼6 3 l  sinðl Þ cosðl Þ
exp ; where we compare the velocity distribution generated in the
n¼1 ln n n n t
simulation with a theoretical normal distribution. For this
ð5Þ
analysis we considered the velocities of the protein carbon
which can be further simplified by retaining only the first term atoms. The linearity of the Q–Q plot indicates that our
in the series,39 which represents a good approximation for simulation data are normally distributed. Application of the
t > 0.2 t, Shapiro–Wilk method (using 4000 samples)40,41 did not show
evidence of deviations from normality.
 2
TðtÞ  Tf 1 ðsinðl1 Þ  l1 cosðl1 ÞÞ2 l1 t
¼6 3 exp ; ð6Þ Computer simulations of the proteins
Ti  Tf l1 l1  sinðl1 Þ cosðl1 Þ t
We have investigated three protein structures, myoglobin
with
(1MBS entry in the Protein Data Bank),42 green fluorescence
GR protein (1QXT),43 and two different states of the Ca2+-
1  l1 cotðl1 Þ ¼ ¼ Bi: ð7Þ ATPase (1SU44 and 1KJU44). The Ca2+-ATPase enzymes
k
are P-type ion pumps that feature in active transport of
Eqn (6) is more convenient for the analysis of the simulation
Ca2+ across the cellular membrane. The catalytic cycle in-
results and is the one we have employed in this work. From the
volves several conformations. 1SU4 is a conformation where
fitting of the simulation results to eqn (6) we can obtain both
Ca2+ is present (E1 state), whereas 1KJU is a conformation
the l1 parameter appearing in eqn (7) and t. The fitting is
corresponding to the Ca2+-free (E2) state. It has been recently
performed as follows. We rewrite eqn (6) as:
found that ATPases release heat when they are transporting
TðtÞ  Tf Ca2+,8 hence knowledge of the thermal transport properties is
ln ¼ A  t=B ð8Þ needed to characterize the heat transfer in these enzymes.
Ti  Tf
Snapshots of the proteins are given in Fig. 2. It can be seen
where A includes the sin and cos terms and B = t/l21. The that Ca2+-ATPases undergo dramatic conformational
fitting of the temperature relaxation data to eqn (8) provides changes during the enzymatic cycle.

1612 | Phys. Chem. Chem. Phys., 2010, 12, 1610–1617 This journal is 
c the Owner Societies 2010
 View Online

This droplet was large enough to solvate the different proteins

without affecting their conformation. Water was modeled
using the simple point charge-extended model (SPC/E).46 All
bonds in the proteins and water were kept rigid using the
LINCS algorithm.47 The van der Waals interactions were cut
off at 1.7 nm. The Coulombic interactions were truncated
at 1.5 nm using a switching function at 1.3 nm. The non-
equilibrium simulations of the protein–water system were
performed without periodic boundary conditions, and both
the linear and angular momenta were removed.48 To avoid the
drift of the protein towards the droplet interface, we reset the
center of mass motion of the protein and solvent separately.
We note that the rotational motion of the protein and solvent
Published on 11 January 2010 on http://pubs.rsc.org | doi:10.1039/B918607G

are not coupled in the simulation timescale investigated in this

work. Considering the protein as a sphere moving in a
newtonian fluid we can estimate the relaxation time associated
to the rotation, t = 4pZR3/(kBT). Using the viscosity of
Downloaded by University of Memphis on 08 October 2012

water, E103 Pa s, the protein radius of gyration, E1 nm

and T = 300 K results in relaxation times of the order of 25 ns.
Fig. 1 (Left) Normalized probability distribution of velocities for This is of the same order of magnitude as experimental
myoglobin. (Right) Quantile–quantile plot to test the normality estimates, 45 ns for myoglobin,49 and much larger than the
of the velocity distribution function; simulations (y-axis), normal simulation times, 0.1 ns, employed to investigate the
distribution (x-axis). temperature relaxation process. Depending on the protein size
the simulations involved 104 (for 1MBS and 1QXT) to 105
(for 1SU4 and 1KJU) solvent molecules, and a time step of 2 fs.
Due to the low vapor pressure of water at 300 K, we did not
observe any significant water evaporation during the simula-
tions. Using this set up the root mean square deviation of the
proteins in the water sphere with respect to the crystallo-
graphic structures was found to be of the order of 0.3 nm.
The transient non-equilibrium simulations were performed
as follows. Firstly we equilibrated the whole system at
constant temperature, 300 K. After this initial equilibration
two Berendsen thermostats50 were applied, one to the protein
and the second one to the water droplet. These thermostats
were applied for about 10 ps to allow the temperatures of the
proteins and water to fluctuate around the corresponding
target values Ti and Tf, respectively. After the equilibration
period, the thermostat on the protein was removed and the

Fig. 2 Snapshots of the proteins investigated in this work. Myoglobin

(top-left), GFP (top-right), 1SU4 (bottom-left) and 1KJU (bottom-
right).

The protein interactions were modeled using the GROMOS96

43a2 force field.45 The structure of the proteins was firstly
minimized in vacuo using the steepest descent method. In order
to perform the transient non-equilibrium simulations, the Fig. 3 Simulation set up employed in the transient non-equilibrium
proteins were subsequently immersed in the center of a simulations, showing the protein Ca2+-ATPase immersed in the water
spherical water droplet with a diameter 9–10 nm (see Fig. 3). droplet.

This journal is 
c the Owner Societies 2010 Phys. Chem. Chem. Phys., 2010, 12, 1610–1617 | 1613
 View Online

Table 1 Heat capacity of the proteins investigated in this work in

kJ mol1 K1. Simulations T = 300 K, experiments 298 K

Protein This work Experiments54 Simulations23

Myoglobin 27  9 24.2 17
GFP 47  10 — 24
Ca2+ATP-ase (1SU4) 180  35 — —
Ca2+ATP-ase (1KJU) 180  30 — —

property has been analyzed by Gomez et al.53 using a wide

range of proteins in different conformations. In that work it
was established that the heat capacity of globular proteins in
solution is dominated, about 80% of the total, by intra-
Published on 11 January 2010 on http://pubs.rsc.org | doi:10.1039/B918607G

molecular contributions, including stretching, bending and

Fig. 4 Variation of the radius of gyration with time for three proteins internal rotations, with intermolecular degrees of freedom
investigated in this work. The simulation results correspond to a contributing only E3%. For proteins in solution it was
transient non-equilibrium simulation, where the initial temperature concluded that there is a hydration contribution to the heat
Downloaded by University of Memphis on 08 October 2012

was set to 350 K. The radius of gyration for two of the proteins has capacity which amounts to about 15% of the heat capacity of
been shifted upwards to facilitate the comparison of the results. the proteins in their native states. Table 1 compares our results
with experimental data available in the literature, which were
system was allowed to relax. The temperature of the water bath obtained using calorimetric techniques. We have also included
was maintained constant during the relaxation process using a simulation data from Yu and Leitner,23 who used a computa-
temperature coupling constant of 0.1 ps. We performed several tional approach to estimate the heat capacities of myoglobin
tests to assess the effect of the temperature coupling constant on and GFP. Our results for the heat capacity of myoglobin are in
the protein relaxation process. We found that the results were good agreement with the experimental results of this protein in
insensitive within statistical uncertainty to the value set for the solution. The estimates from Yu and Leitner are slightly below
coupling constant. We performed an additional test using a the experimental result. We note that these simulations were
Nosé-Hoover (NH) thermostat to assess the influence of the performed in anhydrous conditions.
thermostat. We found that both the Berendsen and NH Our results show a significant increase in the heat capacity
reproduced the same relaxation behavior for the temperature. in moving from myoglobin to GFP and to the Ca2+-ATPases,
The relaxation process was investigated for about 100 ps. The both containing a larger number of carbon atoms, E103 and
equipartition principle expression was employed to compute the E5  103, respectively, as compared with myoglobin
temperature of the protein as a function of time. The temperature E8  102. The heat capacity is an extensive property for
was then fitted to eqn (6) using the Levenberg–Marquadt macroscopic systems, but for small systems consisting of
algorithm.51 The results reported below were obtained from an a few atoms the heat capacity is non extensive.55 Hence our
analysis performed over three independent runs. computations of the heat capacity of single proteins can
The computation of the thermal conductance requires provide a test to assess the extensivity of this property in
previous knowledge of the effective radius of the proteins proteins of varying molecular weight. A representation of the
(radius of gyration) and the heat capacity (see eqn (9)). We experimental heat capacity of the protein as a function of
found that the radius of gyration of the protein fluctuates the number of carbon atoms offers an approach to test the
around a well defined average during the 100 ps duration of consistency of the simulated heat capacities reported above
the relaxation process (see Fig. 4). (see Fig. 5). We note that more involved approaches are
The heat capacity of the proteins was estimated using possible. Gomez et al.53 considered the molecular weight and
constant temperature runs of the protein immersed in water the accessible surface of the proteins, in order to fit heat
using in this case periodic boundary conditions. The heat capacities of proteins in the non-native state. Fig. 5 clearly
capacity was calculated from direct numerical differentiation shows that the experimental results follow a linear dependence
of the internal energy of the protein using temperatures with the number of carbon atoms, as would be expected from
around the temperature of interest, T = 300 K. Typically the extensive character of the heat capacity. Our simulation
we used 5 temperatures, 310, 305, 300, 295 and 290 K, to results conform to this linear behavior and in general are in
calculate the derivative of the internal energy. Each simulation good agreement with the experimental measurements. To the
was 1.2 ns long. All simulations were performed using the best of our knowledge there are no experimental results of the
GROMACS simulation package.52 heat capacity of Ca2+-ATPases. We find that the heat capacity
of these proteins is of the order of 180 kJ mol1 K1. This
result agrees well with the heat capacity obtained from the
Results extrapolation of the experimental results for smaller proteins
(cf. Fig. 5). Within the statistical accuracy of our computations
Heat capacity of proteins
the heat capacities of the two Ca2+ATPase conformations,
The heat capacity plays an important role in determining 1SU4 and 1KJU, are the same, indicating the heat capacity is
the energetics associated to protein folding processes. This not sensitive to the structural detail of the proteins, and

1614 | Phys. Chem. Chem. Phys., 2010, 12, 1610–1617 This journal is 
c the Owner Societies 2010
 View Online
Published on 11 January 2010 on http://pubs.rsc.org | doi:10.1039/B918607G

Fig. 7 Temperature relaxation of the four proteins investigated in

Fig. 5 Heat capacity of native proteins in solution as a function of the this work: myoglobin (top-left), GFP (top-right), ATPase-1SU4
Downloaded by University of Memphis on 08 October 2012

number of carbon atoms in the protein. Experiments 298 K, simula- (bottom-left) and ATPase-1KJU (bottom-right). Ti and Tf represent
tions 300 K. The experimental values (full circles) are taken from the initial and solvent temperatures, 350 and 250 K, respectively. The
ref. 53. The labels denote different proteins; bovine pancreatic trypsin dashed line represents the fitting of the temperature relaxation to
inhibitor (BPTI), ubiquitin (UB), rybonuclease T1 (RT1), lysozyme eqn (6).
(LYS), staphylococcal nuclease (SN), myoglobin (MY), bovine
chymotrypsinogen A (BCA), green fluorescence protein (GFP) and
conductivity and the thermal conductance of the protein–
Ca2+-ATP-ases. The squares represent the simulation results obtained
in this work and the line is a linear fitting to the experimental data.
water interface.
Fig. 7 shows the temperature relaxation results along with
the corresponding fitting using eqn (6). Our model (eqn (6)),
suggesting that the number of degrees of freedom for both which takes into account the different thermal conductivities
conformations is very similar. of the solvent and the protein, provides a good representation
of the temperature relaxation. This is particularly clear for the
Thermal conductance and thermal conductivity of proteins larger proteins (ATPase) where the temperature fluctuations
Fig. 6 shows a representative temperature relaxation curve are significantly reduced as compared to the smaller proteins
for one of the proteins investigated in this work (Ca2+- (cf. Fig. 7).
ATPase-1SU4). All the relaxation processes investigated here Table 2 contains numerical data of the heat transport
involved a time span of about 100 ps. We find that the time properties of the proteins for (Ti, Ts) = (350, 250) K. To test
scale for temperature relaxation is tens of picoseconds. the sensitivity of our results to the simulation conditions we
Assuming, as a first approximation, an exponential decay of performed simulations with different initial and solvent
the temperature with time, exp(t/t) we can estimate a temperatures, (400, 300) K and (350, 250) K. We obtained
characteristic relaxation time, t, for the cooling process. This the same results within the statistical accuracy of our method.
is of the order of 10–20 ps for the proteins investigated in this Hence we report in Table 2 the results for (350, 250) K.
work. We have fitted the simulated temperature relaxation The thermal diffusivities obtained in this work are of the
data to eqn (6) to quantify the thermal diffusivity, thermal order expected for biological tissue56 and they are also in the
range of the thermal diffusivities of small alkane droplets in
water 5–12 Å2 ps1. Our results for the smaller proteins,
myoglobin and GFP, E5 Å2 ps1 are similar to the thermal
diffusivity of Rhodospseudomonas viridis, 7 Å2 ps1, obtained
by different authors14 using a transient non-equilibrium
approach similar to ours. Yu and Leitner23 have employed
more recently a different approach, based on the scaling of the

Table 2 Heat transport properties; thermal diffusivity (D) in Å2 ps1,

thermal conductivity (k) in W K1 m1, thermal conductance (G) in
MW K1 m2), radius of gyration (R) in nm, density (r) in kg m3,
and the Kapitza length (lK) in nm, of the proteins investigated in
this work. All results correspond to temperature relaxations from
Ti = 350 K to Tf = 250 K

Protein D k G R lK
Myoglobin 41 0.13  0.02 100  5 1.529  0.003 1.3
Fig. 6 Temperature relaxation of Ca2+-ATPase (1SU4) as a function GFP 31 0.12  0.01 270  20 1.667  0.002 0.4
of time. The dashed line represents the fitting of the temperature ATPase-1SU4 17  2 0.22  0.02 260  70 3.821  0.004 0.9
ATPase-1KJU 16  2 0.23  0.03 210  40 3.69  0.01 1
relaxation to eqn (6).

This journal is 
c the Owner Societies 2010 Phys. Chem. Chem. Phys., 2010, 12, 1610–1617 | 1615
 View Online

energy diffusion coefficient with the vibrational mode frequency advanced in the discussion of the Biot number (cf. eqn (4)).
of the protein, to estimate the thermal diffusivity of myoglobin The Biot number can be rewritten in terms of the Kapitza
and GFP. These authors equilibrated the proteins in water length as, Bi = R/lK. Hence, for Bi E 1 ) lK E R, and both
before computing the transport properties in vacuo. Our the thermal conductance and the thermal conductivity must be
results lie below the ones reported by these authors, 21.1 included in the description of the thermal relaxation.
and 18.7 Å2 ps1 for GFP and myoglobin, respectively,
although they are comparable in magnitude.
Conclusions and final remarks
The thermal conductivities of the different proteins
are comparable. We obtain values of the order of We have performed an investigation of the thermal properties
0.15–0.2 W K1 m1. These values are similar to the ones of three structurally different proteins; myoglobin, GFP and
obtained in small hydrocarbon droplets,30 and again they are Ca2+-ATPase in water. Myoglobin and ATPases control
comparable to the thermal conductivity reported by Yu and different biological processes, oxygen transport, and active
Leitner, 0.27 W K1 m1 using a different simulation transport of Ca2+ across the cellular membrane, respectively.
Published on 11 January 2010 on http://pubs.rsc.org | doi:10.1039/B918607G

approach. Our results show that all the proteins investi- These proteins are dramatically different in structure and size,
gated in this work are worse heat conductors than water 17951 vs. 109706 Da. In order to obtain information on the
(k E 0.6 W K1 m1) and consequently they should develop transport properties of these biomolecules we have performed
larger temperature gradients as compared with bulk water. transient non-equilibrium molecular dynamics simulations,
Downloaded by University of Memphis on 08 October 2012

The analysis of the temperature relaxation using the whereby we model the temperature relaxation of the protein
solution of the heat diffusion equation, (cf. eqn (6)) provides towards the temperature of the surrounding bath. The analysis
a route to estimate the thermal conductance of the protein– of the temperature relaxation with time offers an approach to
water interface. We note that the protein surface is hetero- extract the thermal conductivity and thermal diffusivity
geneous in chemical composition, and irregular in shape, of the protein as well as the thermal conductance of the
consisting of regions of different curvature. Because the protein–water interface. With this purpose we have considered
complex surface topography of the proteins, consisting of a solution of the diffusion equation that incorporates both the
curved and roughly planar areas, our results for the thermal thermal conductivity of the protein and the thermal conductance
conductance have to be interpreted as illustrative of the as fitting parameters.
average or effective thermal conductance of the protein–water We find that the thermal conductivity of the proteins is low,
interface. Our thermal conductances of the protein–water 0.1–0.2 W K1 m1 as compared with the thermal conductivity
interfaces 100–250 MW K1 m2, are similar to the ones of water, 0.6 W K1 m1, hence proteins should be able to
reported by us for hydrocarbon nanodroplets immersed in sustain large thermal gradients across their structure. This may
water, which are higher than the thermal conductances of have implications on biological processes involving molecular
planar alkane–water interfaces, 65 MW K1 m2.37 In motors such as Ca2+-ATPase. Following our estimates of the
particular, our estimate of the thermal conductance of GFP thermal conductivity, it seems feasible that these proteins are
is in excellent agreement with recent results reported for this capable of maintaining temperature gradients under stationary
protein using a different simulation approach, (i.e. stationary conditions. According to the non-equilibrium thermo-
non-equilibrium simulations), and a different protein force- dynamics theory31–33 the development of such temperature
field (i.e. full atom model).28 These thermal conductances gradient could then, due to reciprocity of phenomena, trigger
are of the order of the ones observed in water–hydrophilic uptake of Ca2+ as well as ATP-hydrolysis, pointing to a
interfaces,57 indicating a strong thermal coupling between possible role of the Ca2+-ATPase as a temperature regulator,
water and the protein surface. We have recently suggested as proposed by de Meis and coworkers.58
that the thermal conductances of curved interfaces can be Within the statistical uncertainty of our computations we do
higher than those of planar interfaces.30 This is a physical not find significant differences in the thermal conductivities of
effect that goes beyond the chemical nature of the surface. the three proteins investigated in this work. The thermal
Although it is not possible to disentangle surface topographic diffusivity seems to be more sensitive to the protein nature
(curvature) and chemical (hydrophilic and hydrophobic though, with an increase in this property of about three times
patches) effects easily, we expect that similar curvature effects in going from the smaller proteins (myoglobin and GFP) to
may apply to the protein–water interface, and hence, the the larger ones (Ca2+-ATPase). The thermal conductances of
corresponding thermal conductance may be higher than that the three proteins vary from 100 to 270 MW K1 m2,
of the corresponding planar interface. myoglobin showing the lower thermal conductance of all the
Finally, using our simulation results we have estimated the systems investigated. Our results for the heat capacity show
Kaptiza length, which represents the thickness of a layer with good agreement with available experimental data of proteins
thermal conductivity, k, which has thermal conductance, G, obtained from calorimetric measurements. In addition to the
i.e., the same as the protein–water interface. The typical values proteins listed above we performed simulations of the heat
for hydrophobic interfaces are of the order of 10 nm, as capacity of lysozyme, rybonuclease T1, and bovine chymo-
compared with 1 nm for hydrophilic interfaces.57 We find trypsinogen (BCA). All are in good agreement with the
Kapitza lengths of the order of E1 nm, i.e., of the order of the experiment, with BCA showing the largest deviations from
radius of gyration of the proteins (see Table 2), showing that the experimental data. The heat capacity correlates linearly with
the interfacial conductance is important for defining the the number of carbon atoms of the proteins, as expected for an
thermal transport of the biomolecules. This idea has been extensive property. Using this result and the experimental data of

1616 | Phys. Chem. Chem. Phys., 2010, 12, 1610–1617 This journal is 
c the Owner Societies 2010
 View Online

small proteins, we have estimated the heat capacity of the Ca2+- 18 D. M. Leitner, Phys. Rev. Lett., 2001, 87, 188102.
ATPase. The extrapolated result is in excellent agreement with 19 D. E. Sagnella, J. E. Straub and D. Thirumalai, J. Chem. Phys.,
2000, 113, 7702–7711.
our computational estimate of 180 kJ mol1 K1. To the best of 20 D. E. Sagnella and J. E. Straub, J. Phys. Chem. B, 2001, 105,
our knowledge this is the first computational estimate of the heat 7057–7063.
capacity of this important protein. 21 K. Moritsugu, O. Miyashita and A. Kidera, J. Phys. Chem. B,
Our work, in conjunction with previous investigations, 2003, 107, 3309–3317.
22 X. Yu and D. M. Leitner, J. Phys. Chem. B, 2003, 107, 1698–1707.
offers a consistent image of the heat transport in biomolecules. 23 X. Yu and D. M. Leitner, J. Chem. Phys., 2005, 122, 054902.
The thermal conductance of the water–protein interface 24 N. Ota and D. A. Agard, J. Mol. Biol., 2005, 351, 345–354.
should play a relevant role in controlling the energy relaxation 25 T. Ishikura and T. Yamato, Chem. Phys. Lett., 2006, 432, 533–537.
26 K. Sharp and J. J. Skinner, Proteins: Struct., Funct., Bioinf., 2006,
of the protein. Overall, our results show that the water–protein 65, 347–361.
interface has a high conductance, similar to the one measured 27 M. Takayanagi, H. Okumura and M. Nagaoka, J. Phys. Chem. B,
in water adsorbed at hydrophilic surfaces.57 We suggest that in 2007, 111, 864–869.
addition to the heterogeneous chemical structure of the 28 N. Shenogina, P. Keblinski and S. Garde, J. Chem. Phys., 2008,
129, 155105.
Published on 11 January 2010 on http://pubs.rsc.org | doi:10.1039/B918607G

proteins, with hydrophobic and hydrophilic patches, other 29 M. J. Assael, E. Charitidou, C. A. N. de Castro and
factors such as the surface topography, characterized by highly W. A. Wakeham, Int. J. Thermophys., 1987, 8, 663–670.
curved interfaces, may be contributing to the enhancement of the 30 A. Lervik, F. Bresme and S. Kjelstrup, Soft Matter, 2009, 5,
thermal conductance of these nanoscale interfaces. We believe 2407–2414.
Downloaded by University of Memphis on 08 October 2012

31 S. Kjelstrup, J. M. Rubi and D. Bedeaux, Phys. Chem. Chem.

this information is relevant to construct microscopic models to Phys., 2005, 7, 4009–4018.
describe how the energy is redistributed in proteins and 32 D. Bedeaux and S. Kjelstrup, Phys. Chem. Chem. Phys., 2008, 10,
through the water–protein interface. In the long term, these 7304–7317.
33 S. Kjelstrup, D. Barragan and D. Bedeaux, Biophys. J., 2009, 96,
fundamental studies may be useful to understand from a 4376–4386.
microscopic perspective certain human diseases, e.g., muscular 34 T. Ikeshoji and B. Hafskjold, Mol. Phys., 1994, 81, 251–261.
dystrophy and sickle cell disease, as well as temperature 35 F. Bresme, B. Hafskjold and I. Wold, J. Phys. Chem., 1996, 100,
regulation in muscle tissue, which have been linked to the 1879–1888.
36 F. Bresme, J. Chem. Phys., 2001, 115, 7564–7574.
Ca2+-ATPase protein.58,59 37 H. A. Patel, S. Garde and P. Keblinski, Nano Lett., 2005, 5,
2225–2231.
38 S. Kjelstrup and D. Bedeaux, Non-equilibrium Thermodynamics of
Acknowledgements Heterogeneous Systems 2008, Series on Statistical Mechanics,
World Scientific, Singapore, Singapore, vol. 16, 2008.
We would like to thank the Imperial College High Performance 39 F. P. Incropera, D. P. DeWitt, T. L. Bergman and A. S. Lavine,
Computing Service and NOTUR (The Norwegian metacentre Fundamentals of Heat and Mass Transfer, John Wiley & Sons,
for computational science), for providing computational Hoboken, USA, 6th edn, 2006.
40 S. Shapiro and M. Wilk, Biometrika, 1965, 52, 591–611.
resources. Funding from the Center for Advanced Study, at 41 P. Royston, Statistical Algorithms, 1995, 44, 547–551.
the Norwegian Academy of Sciences and Letters is gratefully 42 H. Scouloudi and E. Baker, J. Mol. Biol., 1978, 126, 637–660.
acknowledged. 43 D. Barondeau, C. Putnam, C. Kassmann, J. Tainer and E. Getzoff,
Proc. Natl. Acad. Sci. U. S. A., 2003, 100, 12111–12116.
44 C. Xu, W. J. Rice, W. He and D. L. Stokes, J. Mol. Biol., 2002,
References 316, 201–211.
45 L. D. Schuler, X. Daura and W. F. van Gunsteren, J. Comput.
1 P. K. Agarwal, J. Am. Chem. Soc., 2005, 127, 15248–15256. Chem., 2001, 22, 1205–1218.
2 J. F. Swain and L. M. Gierasch, Curr. Opin. Struct. Biol., 2006, 16, 46 H. J. C. Berendsen, J. R. Grigera and T. P. Straatsma, J. Phys.
102–108. Chem., 1987, 91, 6269–6271.
3 J. Ross, J. Phys. Chem. B, 2006, 110, 6987–6990. 47 B. Hess, H. Bekker, H. Berendsen and J. Fraaije, J. Comput.
4 C. Toyoshima, M. Nakasako, H. Nomura and H. Ogawa, Nature, Chem., 1997, 18, 1463–1472.
2000, 405, 647–655. 48 S. C. Harvey, R. K.-Z. Tan and T. E. Cheatham-III, J. Comput.
5 S. W. Lockless and R. Ranganathan, Science, 1999, 286, 295–299. Chem., 1998, 19, 726–740.
6 E. R. Henry, W. A. Eaton and R. M. Hochstrasser, Proc. Natl. 49 G. South and E. Grant, Proc. R. Soc. London, Ser. A, 1972, 328,
Acad. Sci. U. S. A., 1986, 83, 8982–8986. 371–387.
7 V. Botan, E. H. G. Backus, R. Pfister, A. Moretto, M. Crisma, 50 H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren,
C. Toniolo, P. H. Nguyen, G. Stock and P. Hamm, Proc. Natl. A. DiNola and J. R. Haak, J. Chem. Phys., 1984, 81, 3684–3690.
Acad. Sci. U. S. A., 2007, 104, 12749–12754. 51 W. H. Press, S. A. Teukolsky, W. T. Vetterling and B. P. Flannery,
8 M. Suzuki, V. Tseeb, K. Oyama and S. Ishiwata, Biophysical The Art of Scientic Computing, 3rd ed., Cambridge University
Journal: Biophysical Letters, 2007, 92, L46–L48. Press, Cambridge, UK, 2007.
9 F. Bresme, A. Lervik, D. Bedeaux and S. Kjelstrup, Phys. Rev. 52 D. van der Spoel, E. Lindahl, B. Hess, G. Groenhof, A. Mark and
Lett., 2008, 101, 020602. H. Berendsen, J. Comput. Chem., 2005, 26, 1701–1718.
10 D. M. Leitner, Annu. Rev. Phys. Chem., 2008, 59, 233–259. 53 J. Gomez, V. J. Hilser, D. Xie and E. Freire, Proteins: Struct.,
11 Y. Mizutani and T. Kitagawa, Science, 1997, 278, 443–446. Funct., Genet., 1995, 22, 404–412.
12 R. H. Austin, A. Xie, L. van der Meer, B. Redlich, 54 P. Privalov and G. Makhatadze, J. Mol. Biol., 1990, 213, 385–391.
P. A. Lindgoard, H. Frauenfelder and D. Fu, Phys. Rev. Lett., 55 T. L. Hill, Thermodynamics of small systems, Dover, Mineola,
2005, 94, 128101. New York, 2002.
13 R. J. D. Miller, Annu. Rev. Phys. Chem., 1991, 42, 581–614. 56 Y. Touloukian, R. Powell, C. Ho and M. Nicolaou, Thermo-
14 M. Tesch and K. Schulten, Chem. Phys. Lett., 1990, 169, 97–102. physical Properties of Matter, Plenum Press, New York, 1973.
15 P. Li and P. M. Champion, Biophys. J., 1994, 66, 430–436. 57 Z. Ge, D. G. Cahill and P. V. Braun, Phys. Rev. Lett., 2006, 96,
16 K. Blumhagen, I. Muegge and E. W. Knapp, Int. J. Quantum 186101.
Chem., 1996, 59, 271–279. 58 L. de Meis, G. Oliveira, A. Arruda, R. Santos, R. da Costa and
17 K. Moritsugu, O. Miyashita and A. Kidera, Phys. Rev. Lett., 2000, M. Benchimol, Life, 2005, 57, 337–345.
85, 3970–3973. 59 E. Carafoli, Physiological Reviews, 1991, 71, 129–153.

This journal is 
c the Owner Societies 2010 Phys. Chem. Chem. Phys., 2010, 12, 1610–1617 | 1617