Generic folding and transition hierarchies for surface adsorption of hydrophobic-polar lattice model proteins

1.2 MB·English

by  Ying Wai Li

Generic folding and transition hierarchies for surface adsorption of HP lattice model
proteins
Ying Wai Li∗
Center for Simulational Physics, University of Georgia, Athens, Georgia 30602, U.S.A.†
Thomas Wu¨st‡
Swiss Federal Research Institute WSL, Zu¨rcherstrasse 111, CH-8903 Birmensdorf, Switzerland
David P. Landau§
Center for Simulational Physics, University of Georgia, Athens, Georgia 30602, U.S.A.
3
1 (Dated: January 16, 2013)
0
The thermodynamic behavior and structural properties of hydrophobic-polar (HP) lattice pro-
2
teins interacting with attractive surfaces are studied by means of Wang-Landau sampling. Three

n benchmark HP sequences (48mer, 67mer, and 103mer) are considered with different types of sur-
a faces, each ofwhich attracting eitherall monomers, only hydrophobic(H)monomers, oronly polar
J (P) monomers, respectively. The diversity of folding behavior in dependence of surface strength

5 is discussed. Analyzing the combined patterns of various structural observables, such as e.g., the
1 derivatives of the numbersof surface contacts, together with the specific heat, we are able to iden-
tify generic categories of folding and transition hierarchies. We also infer a connection between

] thesetransition categories andtherelativesurface strengths,i.e., theratio ofthesurface attractive
t
f strength to the inter-chain attraction among H monomers. The validity of our proposed classifi-
o
cation scheme is reinforced by the analysis of additional benchmark sequences. We thus believe
s
that the folding hierarchies and identification scheme are generic for HP proteins interacting with
.
t attractive surfaces, regardless of chain length, sequence,or surface attraction.
a
m
PACSnumbers: 05.10.-a,82.35.Gh,87.15.ak,87.15.Cc
-
d
n
I. INTRODUCTION they occur remain puzzles to be solved [1, 13]. Enor-
o
mous efforts have been dedicated to unveil the mysteries
c
[ Protein adsorption on solid surfaces has attracted in protein folding and adsorption mainly by experimen-

strong research interest recently for its numerous appli- tal approaches [14, 15]. Nevertheless, the fact that only
1
v cations in nanotechnology and biomaterials [1, 2]. The the “final product” can be obtained and studied in an
2 study of protein functions in experiments often involves experimentmakesforslowprogressinunderstandingthe
6 the immobilization of proteins [3, 4]. Adhesion of pro- dynamics and folding processes. From another point of
4 teins on solid substrates including metals, semiconduc- view, the diversity of possible protein sequences and so-
3 tors, carbon or silica etc., enables the synthesis of new phisticated interactions among amino acids also compli-
1. materials for biosensors or electronic devices [5–8]. It is cate theoretical structural prediction in protein folding -
0 also the key to reveal the principles of many biological not to mention when the protein interacts with solvent
3 processes and causes of diseases, e.g., when integrating moleculesorasubstratewhereanextralevelofcomplex-
1 implanted materials with body tissues [9, 10]. Under- ity enters.
:
v standing how the functions and conformations of a pro-
With the advances in computer power,numericalsim-
i tein are affected by adsorption and desorption is an im-
X ulation has become a promising way to shed more light
portant topic in protein drug delivery [11, 12].
r on the general problem. The study of simplified, coarse-
a It is known that configurational changes of protein
grainedproteinmodels in conjunction with Monte Carlo
molecules upon surface adsorption depend on both the
simulations, being able to efficiently explore large con-
protein’s properties (e.g. sequence, size, thermodynamic
formational phase spaces [16, 17], has been a particu-
stability, etc.) and the surface properties (e.g. mate-
larlysuccessfulapproachinunderstandingthe principles
rials, polarity, surface roughness, etc.); but how large
behind protein folding and adsorption. Among these
these changes are and where in the protein molecules
minimalist models, the hydrophobic-polar (HP) model
introduced by Dill et al. [18, 19] has been a particu-
larlyactive,yetdifficult,researchsubjectinrecentyears.
The interest in the HP model is twofold: (i) Despite its
∗ [email protected]
known limitations [20, 21], the model captures some of
† Present address: National Center for Computational Sciences,
the most important qualitative features which drive the
Oak Ridge National Laboratory, Oak Ridge, Tennessee 37831,
U.S.A. folding of proteins and characterize their native states.
‡ [email protected] It thus laid a basis to systematically study many prob-
§ [email protected] lems in protein folding by means of computer simula-
2
tions. (ii)The groundstate searchforanHPsequence is itcanonlybehydrophobic(H)orpolar(P).Aproteinis
anNP-complete problem[22, 23], and the roughfree en- thusrepresentedbyaheteropolymerwhichconsistsofN
ergy landscapes cause traditional Monte Carlo methods, connectedmonomersoftype HorP.Anattractiveinter-
e.g. Metropolissampling [24],to failinthe lowtempera- action exists only between a pair of non-bonded nearest
tureregime[25]. Therefore,theHPmodelisawell-suited neighboring H monomers. This attraction is denoted by
testinggroundforanumberofemergingnumericalmeth- ε in our discussion, and the magnitude indicates the
HH
ods[26–34]. Inaddition,the thermodynamicbehaviorof abilityofthe Hmonomersto pullthemselvestogetheras
different HP sequences can vary noticeably, even for the determined by the insolubility of the protein in an aque-
same chain length [32]. Consequently, finite-size scal- ous environment. In other words, the solvent quality is
ing cannot be applied to a systematic study of the effect intrinsically considered by the model. Other factors like
of system size, unlike many other models in statistical charges and acidity of amino acids that also govern pro-
physics. Thus, our understanding of the general behav- tein folding are not handled in this scope.
ior of the model is still incomplete and the simulation Inadditiontotheinternalinteractionswithinthepoly-
of the HP model remains a challenging computational mer, the binding of a protein with an attractive sub-
problem in statistical physics. strate contributes to the total energy. We have consid-
Various approaches have been undertaken to better ered three types of surface fields in view of the setting
understand the energy landscapes, thermodynamics and of the HP model: (i) a surface interacts only with H
conformational transitions of HP proteins adsorbing on monomers (a hydrophobic surface) with strength ε ,
SH
an attracting surface [35–40]. In this work, our intent is (ii) a surface interacts only with P monomers (a polar
toidentifygenericthermodynamicandstructuralbehav- surface) with strength ε , and (iii) a surface interacts
SP
ior of proteinadsorptionfrom a macroscopicperspective with both H and P monomers with equal strength, i.e.,
using the HP model. We have adopted Wang-Landau ε =ε 6=0. The energy function of an HP sequence
SH SP
sampling [41–43] to obtain the density of states of a sys- interacting with a surface takes the general form:
tem,fromwhichsubsequentthermodynamicscanbecal-
culated. Tothebestofourknowledge,itisthefirstcom- E =−εHHnHH −εSHnSH −εSPnSP, (1)
prehensive analysis of structural transitions for protein
where n denotes the number of interacting pairs be-
adsorption that integrates results from multiple HP se- HH
tween H monomers, n the number of surface contacts
quences. Previous work by Bachmann and Janke [39], SH
with H monomers and n the number of surface con-
Swetnam and Allen [40] or Radhakrishna et al. [44] SP
tacts with the P monomers. Besides contributing to the
only studied the conformational pseudo-phases based on
system’senergy,thesethreequantitiesarealsouseful“or-
individual benchmark HP sequences. Comparing the
der” parameters that give quantitative measures of the
thermodynamic and structural properties of multiple se-
structure of a conformation [45]. The negative signs in
quences allowedus to identify categories of generic tran-
front of each term indicate that it is energetically favor-
sition behavior and to draw a correspondence between
ablewhenthemonomersinteractorcomeincontactwith
thesecategoriesandtherelativeinteractionstrengthsin-
the surface.
volved. This finding is fundamental as it implies that
Our simulations were performed on a three-
different HP sequences share certain general, qualita-
dimensional simple cubic lattice with the attractive
tive, characteristics in structural transformations when
surface being the xy-plane placed at z = 0. A non-
brought near to an adsorbing substrate.
interacting steric wall is placed at z = N +1 to confine
Thearticleisorganizedasfollows: SectionIIdescribes
the polymer in a way that it can contact both walls
the model employed, Section III explains the sampling
with its ends only when it is a vertical straight chain.
method and details, Section IV presents the simulation
Periodic boundary conditions are imposed in the x and
results, Section V discusses the analysis and interpreta-
y directions.
tion of the results, and Section VI gives a summary and
outlook of this work.
III. METHOD
II. MODEL
A. Calculation of thermodynamic quantities
Thecommonlyknown22aminoacidsfoundinproteins
The partitionfunction, Z,ataparticulartemperature
can roughly be classified as either hydrophobic or polar
T can be expressed in terms of the energy density of
depending on the nature of their side chains. The ten-
states g(E):
dency of the non-polar residues to stay away from water
molecules has been identified as the key driving “force” Z = g(E)e−E/kBT, (2)
in forming tertiary structures. The hydrophobic-polar
XE
(HP) model [18, 19] is a coarse-grainedlattice model for
proteins that captures this hydrophobic effect. In this where E is the energy of the system as defined by the
model,anaminoacidistreatedasasinglemonomerand energy function, k is the Boltzmann constant, and the
B
3
sum runs over all the energies that the system can take. whichcanbe performedfromstate A; n andn are
B→A B
Since g(E) does not depend on T, one may calculate defined likewise. Because of reversibility of pull moves,
Z at any temperature with a single estimate of g(E). n = n and the two terms cancel out in Eq.
A→B B→A
All the thermodynamic quantities then follow from the (5). When a pull move is chosen to generate a new con-
knowledge of Z. For example, the average energy hEi figuration, a list of possible moves from state A is first
and the heat capacity C are calculated as: constructed to obtain n . A move is then selected ran-
V A
domly from the list, generating state B. n is counted
B
1
hEi= Eg(E)e−E/kBT, (3) and P(A → B) can then be calculated. This procedure
Z XE is computationally expensive (it slows down the simula-
tion by approximately an order of magnitude); however,
it secures the reliability of our results by taking detailed
E2 −hEi2 balance into careful consideration.
C = . (4)
V (cid:10) k(cid:11) T2 We adopted RANLUX as the random number gener-
B
ator (gsl rng ranlxd2), which is recommended by the
The specific heat is then defined as C /N accordingly.
V GNU Scientific Library (GSL) for its “reliable source of
uncorrelated numbers” and “strongest proof of random-
ness”[51]. Itusesalagged-Fibonacci-with-skippingalgo-
B. Wang-Landau sampling and trial moves
rithm[52]withdoubleprecisionoutputandalongperiod
171
of about 10 .
Wang-Landau (WL) sampling is a powerful, iterative
The use of Wang-Landau sampling together with in-
algorithm to estimate the density of states, g(E), with
ventivetrialmovesisessential foruncoveringthe subtle,
high accuracy. Details of the algorithm are described in
lowtemperaturethermodynamicsofthesesystems. This
Refs.[41–43]. Inoursimulations,ratherstringentparam-
is demonstrated in Fig. 1 which shows the specific heat
eters were chosen in order to obtain accurate estimates
of a 36mer (P3H2P2H2P5H7P2H2P4H2P2HP2) interact-
for g(E): We used a flatness criterion p = 0.8 for the
ing with a very weak attractive surface (ε = ε =
SH SP
48mer and p=0.6 for the 67mer and 103mer. The final
1,ε =12). Asseeninthefigure,allthesurfacerelated
modification factor was set to ln(f ) = 10−8 in all HH
final transitions below k T/ε ≈ 0.3 are completely inac-
B HH
cases.
cessiblebyMetropolissampling,evenwithverylongruns
It has been found that two types of trial moves, pull (108 trialmoves!) incorporatingpullmovesandbondre-
moves [46] and bond-rebridging moves [47], work partic-
bridging moves! In contrast, Wang-Landau sampling is
ularly well together with WL sampling in search of the
even able to uncover the shoulder at k T/ε ≈ 0.016
B HH
global energy minimum conformations and the determi-
shown in the inset (caused by a rapid change in configu-
nation of the density of states for lattice polymers [48–
rational entropy of the excitation from the ground state
50]. Theabilityofreachinglowenergystatesallowsfora
more thorough survey of conformational space, yielding
a higher resolution of g(E) and thus more precise ther-
Cv / N
modynamic quantities especially in the low temperature
1.2
regime. This is of particular importance for longer chain
0.2
lengths with more complex energy landscapes [50].
The two trialmoves are calledwith different probabil-
ities. Bond-rebridging moves transform a polymer from 0.1
one compact state to another without uncoiling, making
it more efficient than pull moves in dealing with densely
packed polymers. However, since the energy difference 0.00 0.02 0.04
resultingfromabond-rebridgingmoveisrelativelylarge,
0.4
its acceptance rate is rather low. This drawback is com-
Wang-Landau
pensated for with a higher calling ratio. In our simula-
Metropolis
tions, we used move fractions of 80% and 20% for bond-
rebridging moves and pull moves, respectively.
Inordertofulfilldetailedbalancewhenemployingpull 0.0
0.0 0.2 0.4 0.6
moves,anextrafactoris addedtothe acceptanceproba- k T / ε
B HH
bility of moving from state A (with energy E ) to state
A
B (with energy EB): FIG.1. (Coloronline)Maingraph: Specificheatofa36mer
interactingwithaveryweakattractivesurfaceasobtainedby
g(EA)nB→A/nB Wang-Landau and Metropolis sampling. Where not shown
P(A→B)=min 1, , (5)
(cid:18) g(E )n /n (cid:19) statistical errors are smaller than the data points. Compa-
B A→B A
rable computational effort was used for each method. Inset:
where n is the number of pull moves that trans- Resultsatextremelylowtemperature. Typicalconfigurations
A→B
form A to B; n is the number of possible pull moves aboveand below the shoulderare shown.
A
4
to the first few excited states; see Ref. [53] for details). and
We, therefore, stress that sophisticated Monte Carlo al-
1
gorithms are crucial in simulating systems with delicate hQi= Qg(E,Q)e−E/kBT. (10)
Z
low temperature behavior. Q EX,Q
Recently, improvements have been proposed to speed
up and ensure the convergenceof WL sampling in simu-
lating lattice polymers or proteins [40, 44, 50, 54]. IV. SIMULATION RESULTS
We have studied three benchmark HP sequences
C. Calculation of thermodynamics of structural (48mer,67mer,and103mer)interactingwiththreetypes
observables ofsurfaces(seeabove). The48mer(PHPHP4HPHPHP2-
HPH6P2H3PHP2HPH2P2HPH3P4H) is seq. 9 among
Structural quantities are essential in understanding the ten “Harvard sequences” designed originally for
non-energeticproperties of the system. They providein- algorithm testing purpose [57], where the number
formation on the structures and packing of the polymer. of H and P monomers are exactly the same in each
Apartfromthequantitiesn ,n andn introduced sequence, i.e., 24 H monomers and 24 P monomers
HH SH SP
in Section II, structural quantities that are often of in- respectively. The one chosen here for our simulation
terest include the radius of gyration, has the minimum ground state degeneracy in 3D
free space [32, 57]. The 67mer (PHPH2PH2PHPP-
N H3P3HPH2PH2PHP2H3P3HPH2PH2PHP2H3P3HP-
1
R =v (~r −~r )2, (6) H2PH2PHP2H3P) was first introduced to resemble
g uuN Xi=1 i cm α/β-barrel in real proteins [58], while the 103mer
t (P2H2P5H2P2H2PHP2HP7HP3H2PH2P6HP2HPHP2H-
and the end-to-end distance, P5H3P4H2PH2P5H2P4H4PHP8H5P2HP2) was proposed
to model cytochrome c [59].
Ree =|~rN −~r1|, (7)
where ~r in Eq. (6) is the center of mass of the config- A. Ground states and limiting behavior
cm
uration;~r is the position of monomer i.
i
To obtain the thermodynamics of a structural observ- TableIreportsthelowestenergiesfoundforthesethree
ableQ,weestimatedthe jointdensityofstates, g(E,Q), sequences during the estimation of g(E) with different
by multicanonical sampling [55, 56]. Although it is types of attractive surfaces. To understand the “asymp-
feasible to sample g(E,Q) all over again using a two- totic” folding behavior inthe limit ofa surface with infi-
dimensional random walk in WL sampling if only one niteattractivestrength,wesimulatedthe casewherethe
structural quantity is required, it becomes impossible if HP chain was confined to a two-dimensional free space.
severalofthem areofinterest. Amoreefficientwayis to This is equivalent to restricting all monomers of the HP
make use of the prior knowledge of g(E) and perform a chain on the surface, giving a 2D ground state. Another
multicanonical sampling. In this process, trial states are limiting case is when the surface is totally absent in a
accepted or rejected according to 1/g(E), where g(E) is three-dimensional space. For more details on these two
the one obtained previously by the one-dimensional WL cases, see [60].
samplingandisheldfixedthroughoutthe wholeproduc- These two limiting cases are useful in visualizing up-
tion procedure. As a new trial state is accepted, any perandlowerboundsforthermodynamicobservablesand
desired structural quantity Q would be calculated and they serve as an aid to understand the details of folding
accumulated in a two-dimensional histogram, H(E,Q). behavior. A demonstrative example is the averaged ra-
The simulation is brought to an end when a sufficiently dius of gyration per monomer, hR i/N, shown in Fig.
g
largenumberofMonteCarlostepshavebeenperformed. 2 for the 48mer interacting with a surface attracting all
The joint density of states, g(E,Q), is then obtained by monomers (surfaces A1, A2 and A1/2). The radii of gy-
reweighing H(E,Q): ration of the two limiting cases are plotted on the same
figure. Drawing a simple connection to the self-avoiding
g(E,Q)=g(E)H(E,Q). (8) random walk on square and cubic lattices, it is obvious
that hR i is largestwhen allthe monomers are forcedto
g
Assuch,wecanobtainasmanyg(E,Q)’sforvariousQ’s sit on the surface to form planar structures, yielding the
as desired in a single production run. upper bound for hR i. Correspondingly,hR iis smallest
g g
The partition function, Z , for observable Q and its when the HP chain is allowed to fold freely in a three-
Q
expectation value can then be obtained as dimensionalspacetoform3Dstructures,givingthelower
bound of hR i. Generally speaking, when the HP chain
g
Z = g(E,Q)e−E/kBT, (9) is placed near a surface of finite attractive strength, it
Q
EX,Q remainsasanextendedcoilathightemperatureasifthe
5
TABLE I. Lowest energies found for the 48mer, 67mer, and 0.8 Cv
103merinteractingwithdifferentattractivesurfacetypesand
N
strengths,whichareabbreviatedinthesurfacelabels(A,Hor
P standforthesurfacetypes,thenumbersstandfortheratio 0.6
between εSH or εSP and εHH). Classification of transition
categoriesaredenotedbytheRomannumbersinparentheses.
0.4
surface εHH εSH εSP 48mer 67mer 103mer
Freespace without surface: 0.2
2D 1 / / −21 −29 −32
3D 1 / / −34 −56 −58
0.0
Surfaces attract all monomers: <R >
g
A1 1 1 1 −69 (I) −96 (I) −135 (I)
0.12 N
A2 1 2 2 −117 (I) −163 (I) /
A1/2 2 1 1 −93 (II) −132 (II) −167 (I/II)
0.10
Surfaces attract only H monomers:
H1 1 1 0 −49 (II) −72 (II) −80 (II) 0.08
H2 1 2 0 −73 (I) −108 (I) /
2D free space
H1/2 2 1 0 −79 (III) −118 (II) −128 (III) surface A2: ε = 1, ε = ε = 2
0.06 HH SH SP
surface A1: ε = 1, ε = ε = 1
HH SH SP
Surfaces attract only P monomers: surface A½: ε = 2, ε = ε = 1
HH SH SP
P1 1 0 1 −48 (II) −69 (II) −100 (I/II) 0.04 3D free space
P2 1 0 2 −71 (I) −91 (I) /
P1/2 2 0 1 −79 (III) −123 (II) −150 (II) 0.0 1.0 2.0 3.0 4.0 5.0
k T / ε
B HH
surface is absent. The radii of gyrationfor all cases thus FIG. 2. (Color online) Upper panel: Specific heat CV/N
coincide with the 3D, surface-free one. as a function of the effective temperature kBT/εHH. Lower
As the temperature decreases, the HP chain interact- panel: Averaged radii of gyration per monomer, hRgi/N,
ing with astrongerattractivesurface (A2)starts the ad- as a function of kBT/εHH, for the 48mer interacting with a
surface which attracts all monomers with different strengths
sorption process the earliest at k T/ε ≈ 5.0 as its
B HH (surfacesA1,A2andA1/2). NotethatkBT isscaledwiththe
hR i “departs” from the lower bound and begins to ap-
g internal attraction strength, εHH, so as to compare different
proach the upper bound. Such an adsorption “transi-
systemsin thesame energy scale. Inthis manner,any differ-
tion” is clearly signaled by the peak in CV centered at enceinquantitiescomessolelyfromthesurfacestrengthsεSH
kBT/εHH ≈ 4.25. At kBT/εHH ≈ 1.0, hRgi merges and εSP. Errors smaller than thedata pointsare not shown.
with the upper bound signifying a complete adsorption
ofallmonomers. Theformationofahydrophobiccorein
which the number of inter-chain H-H interactions, n ,
HH
is maximized, then takes place entirely on the surface in
this case until the ground state is reached at zero tem-
perature. This process in the low temperature regime
is identical to the one in two-dimensional free space, as
indicatedbythecompleteagreementintheradiiofgyra-
tion and the coincidence of the peak at k T/ε ≈ 0.5
B HH
observed in CV. WhilethetwopeaksinthespecificheatofsurfaceA1/2
The thermodynamics for surface A1 is qualitatively tend to give the impression that it has the same quali-
similar to that of surface A2 except that it requires a tative folding behavior as the previous cases, the shape
lower adsorption temperature. Since the radii of gyra- ofthe radiusof gyrationclearlydistinguishes it fromthe
tion for both surface types end up with the same value others, apart from showing that the ground state is now
as the upper bound at T = 0, one may expect that the three-dimensional. This is the first clue that the specific
ground state conformations for both systems are two- heat alone does not provide a complete picture of struc-
dimensional. This has been confirmed by the number of tural transition behavior. Indeed, the transition hierar-
surface contacts (n = n = 24, meaning the entire chyofthe48merinteractingwithsurfaceA1/2isdifferent
SH SP
chain is in contact with the surface) and the number of from surfaces A1 and A2, which can only be verified by
H-H interactions (n = 21, which is the same as the examining other structuralparameters as we shall see in
HH
ground state of the 2D limiting case). the following section.
6
0.8 0.12
B. Identification of transition categories by
canonical analysis of specific heat and structural
quantities <R g > / N
0.6 0.10
Three major transitioncategories were identified from
0.4 0.08
the 24 systems presented in Table I by considering the
combined patterns of the specific heat C /N and the
V
average radius of gyration hRgi/N. 0.2 Cv / N 0.06
CategoryI:C showstwopeaks,abumpbetweenthe
V
peaksmightbepossible,hR ishowsamaximumbetween
g
these two peaks. 0.0 0.04
Category II: C shows two peaks, hR i decreases
V g
upon cooling. In the very low temperature regime, it
might rise back up a little to form a minimum when the
temperature approaches zero.
Category III: C shows only one peak with possible
V
shoulders, hR i decreases upon cooling. 0.0
g
ThedifferentcombinationsofC andhR iinthetran-
V g
sition categories are caused by the different order of oc-
currence in folding processes, which are revealed by fur-
-0.2
ther investigation of proper structural parameters and 1d<n >
SH
their derivatives. Quantities which are particularly in- N dT
formative for our systems include the derivatives of the
average number of H-H interactions, dhn i/dT, and -0.4
HH
thoseofthenumbersofsurfacecontacts,dhnSHi/dT and 1d<nH H >
dhn i/dT. Apeakindhn i/dT signalstheconstruc- N dT
SP HH
tion of H-H interactions to form a hydrophobic core (H-
-0.6
core formation). Peaks in dhn i/dT and dhn i/dT 0.0 1.0 2.0 3.0 4.0 5.0 6.0
SH SP
provide information about the formation of surface con- k T / ε
B HH
tacts, which is associatedwith the adsorption process as
well as the “flattening” of structure due to surface at-
traction. FIG. 3. (Color online) Thermodynamics of the 67mer in-
We observe that hReei behaves quite similarly as hRgi teracting with surface A2 (εHH = 1,εSH = εSP = 2), which
butislessreliableatlowtemperaturewheremainlycom- showsatypicalcategoryItransition. Upperpanel: Specific
pact structures are found. For this reason our analysis heat,CV/N,andtheaverageradiusofgyrationpermonomer,
relies on hR i rather than on hR i. hRgi/N,asafunctionoftheeffectivetemperaturekBT/εHH.
g ee The horizontal arrows beside the labels indicate the axes to
which the quantities refer. Middle panel: Typical config-
urations at different temperatures. Lower panel: Deriva-
1. Folding behavior with a strong attractive surface: tives of the average numbers of H-H contacts per monomer,
Category I (1/N)dhnHHi/dT,andthatoftheaveragenumberofsurface
contactsofHmonomerspermonomer,(1/N)dhnSHi/dT,as
Figure 3 shows a typical transition pattern of cate- a function of kBT/εHH. Errors smaller than the data points
are not shown.
gory I demonstrated by the 67mer with surface A2, for
which ε = 1,ε = ε = 2. It is characterized
HH SH SP
by two pronounced peaks in C , with hR i attaining
V g
its maximum between them as seen in the upper panel A closer look at C in Fig. 3 shows a weak bump be-
V
of the figure. The nature of transitions to which the tweenk T/ε ≈1.0and3.5. Thesamephenomenonis
B HH
two peaks in C correspond can be identified by com- alsoobservedindhn i/dT (anddhn i/dT),suggest-
V SH SP
paring them with dhn i/dT and dhn i/dT in the ing that a subtle process related to the interaction with
HH SH
lower panel. Since the surface attracts both types of thesurfaceistakingplaceinthistemperaturerange. Re-
monomers equally, dhn i/dT shows similar behavior calling our previous discussion on the radius of gyration
SP
as dhn i/dT and thus is not shown in the figure. The hR iinSectionIVA,thisisaregionwheretheHPchain
SH g
C peak at k T/ε ≈ 4.4 represents the desorption- keeps forming contacts with the surface until it adsorbs
V B HH
adsorption transition where dhn i/dT peaks at the completely on the surface. We call this process a “flat-
SH
sametemperature. TheC peakatk T/ε ≈0.4rep- tening” of the structure. When the surface attraction is
V B HH
resents the H-core formation as dhn i/dT also shows sufficiently strong, it occurs right after the chain is ad-
HH
a peak at that position. hR i decreases most rapidly sorbed to the surface but before the H-core formation.
g
during this latter process when temperature is lowered. The top graph in the leftmost column in Table II shows
7
a similar case for the 48mer, ε =1,ε =ε =2. 1.2 0.08
HH SH SP
However,therearecaseswherethis “flattening”bump
1.0 0.07
is not observedin C , as seen from the other two exam-
V Cv / N <R g > / N
ples for the 67mer and 103mer with different surface at- 0.8 0.06
tractions in Table II. The flattening process might have
0.6 0.05
been “integrated” within adsorption, or it simply does
notcauseenoughenergyfluctuationstogiveavisiblesig-
0.4 0.04
nalinC . Inthelattercase,signalscanbefoundinother
V
structural quantities like dhn i/dT or dhn i/dT. 0.2 0.03
SH SP
0.0 0.02
2. Folding behavior with a moderately attractive surface:
Category II
Figure 4 shows the thermodynamics for the 103mer
with surface P1/2, a typical case in category II. Simi-
lar to category I, systems in category II also show two
pronounced peaks in C and identification of structural
V
transitions depends on the derivatives of hnHHi, hnSHi 0.2 1 d<n >
SH
and hnSPi. The peak at kBT/εHH ≈ 0.85 represents N dT
the desorption-adsorption transition as identified by the 0.0
peaks in dhn i/dT and dhn i/dT. Another peak at
SH SP
kBT/εHH ≈ 0.42 indicates the H-core formation as sig- -0.2
naled by a peak in dhn i/dT.
HH
The feature that differentiates category II from cate- -0.4
goryIistheabsenceofamaximumforhRgibetweenthe 1 d<nH H > 1 d<nS P >
twopeaksinC . Itdecreasesuponcoolinguntilthevery -0.6 N dT N dT
V
low temperature regime. The difference arises from the
0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4
fact that the flattening of structures occurs at a lower k T / ε
B HH
temperature thanthe H-coreformationinthe vicinity of
a less attractivesurface,givingrise to anothertransition
hierarchy than category I. Two possibilities for hR i are
g FIG. 4. (Color online) Thermodynamics of the 103mer in-
then observed when the temperature is further lowered:
teracting with surface P1/2 (εHH = 2,εSH = 0,εSP = 1),
(a) it keeps descending as in Fig. 4; (b) it rises back up
which shows a typical category II transition. Upper panel:
until T = 0, forming a minimum below the H-core for- Specific heat, CV/N, and the average radius of gyration per
mation temperature as the 48mer does in Table II (top monomer, hRgi/N, as a function of the effective tempera-
graph, middle column for category II). turekBT/εHH. Thehorizontal arrows besidethelabels indi-
Interesting observations at low temperature are cate the axes to which the quantities refer. Middle panel:
revealed by the thermodynamics of dhn i/dT, Typical configurations at different temperatures. Lower
HH
panel: Derivatives of the average numbers of H-H contacts
dhn i/dT anddhn i/dT asshowninthelowerpanel
SH SP
of Fig. 4. During H-core formation at k T/ε ≈ 0.42 per monomer, (1/N)dhnHHi/dT, and those of the numbers
where dhn i/dT peaks at, troughs aBre obHsHerved in ofsurfacecontacts,(1/N)dhnSHi/dT and(1/N)dhnSPi/dT,
HH as a function of kBT/εHH, respectively. Errors smaller than
dhn i/dT anddhn i/dT. Thisisaprocessof“thick-
SH SP thedata points are not shown.
ening” during which some of the surface attachments
have to be broken to facilitate the construction of H-H
interactions.
When the temperature is further lowered to
k T/ε ≈ 0.25, a subtle shoulder could barely be
B HH
seen in C and hR i stays still on cooling; dhn i/dT,
V g SP
however, shows a clear peak. This suggests that surface
contacts for the P monomers are established, demon-
strating the flattening effect. Eventually the structures
with minimal possible energy are attained but they Inmanyotherexamplesofthistransitioncategory(see
no longer span as many layers vertically as at higher TableII),C onlyhastwomajortransitionpeaks,some-
V
temperature. Thesestructuresarenotcompletelyplanar times with a subtle shoulder or a spike merged into the
asincategoryI,asformingsurfacecontactsisnotalways peaks as a result of a combination of various events. In-
more energetically favorable than forming hydrophobic dividual investigation of structural measures is thus es-
H-H interactions. sential to segregate different structural changes.
8
1.25
one may identify adsorption at k T/ε ≈ 0.63 and
0.08 B HH
flattening atk T/ε ≈0.24respectively. dhn i/dT
B HH HH
1.00
0.07 manifests a wide peak across the adsorptionand flatten-
Cv / N <R g > / N ing temperatures, which suggests that the hydrophobic
0.75 0.06 core is formed roughly in the temperature range
k T/ε ≈0.25−0.87. Insteadofproducingindividual
B HH
0.50 0.05 peaks in C , the signals of adsorptionandflattening are
V
“bridged” and smoothed out by the H-core formation,
0.04
0.25 giving only a peak with a shoulder in C , the shape
V
0.03 and location of the latter being often subject to some
0.00 variability. Therefore, it is necessary to rely again on
structural quantities to separate signals for the various
transitions as illustrated in the previous categories.
For veryweak attractivesurfaces,adsorptionand flat-
tening occur at even lower temperatures. They become
distinguishable from the H-core formation which takes
placeatahighertemperature,formingtwooreventhree
distinct peaks in C . We generallyclassify systems with
V
0.0 1d<n >
SH this transition hierarchy as category IV. For a detailed
N dT
discussionof an example from this category,see [53, 61].
-0.1
V. DISCUSSION
-0.2 1d<n > 1d<n >
SP HH
N dT N dT
A. Remarks on the structural measures and
-0.3
categories
-0.4
Ourresultsdemonstratethatacomprehensiveanalysis
0.0 0.5 1.0 1.5
k T / of the specific heat (CV) combined with a set of appro-
B HH priate structural quantities is essential to shed light on
recognizing structural transformations, especially those
FIG. 5. (Color online) Thermodynamics of the 48mer inter-
subtle ones for which C alone provides insufficient in-
acting with surface P1/2 (εHH =2,εSH =0,εSP =1), which formation. We also noteVthat in identifying phase tran-
shows a typical category III transition. See Fig. 4 for figure
sitions, the peaks observed in structural quantities and
explanations.
those in C might be slightly off. One possible expla-
V
nation could be finite size effects [62]. Nevertheless, this
does not affect our identification scheme much as these
3. Folding behavior with a weak attractive surface:
shifts are sufficiently small comparedto the difference in
Category III and beyond
temperature scales required to clearly identify distinct
categories of transition patterns.
When the surface attractive strength further reduces, Table II further summarizes and compares the ther-
the adsorption and flattening temperatures decrease ac- modynamic features for the three categoriesdiscussedin
cordingly. CategoryIII is identified whenthe adsorption theprevioussection. Wehaveintentionallychosenexam-
transition coincides with H-core formation, giving a sin- plesofvariouscombinationsofsurfacetypes(i.e. surface
gle peak in C associated with a shoulder in some cases strengths) and chain lengths for each category to illus-
V
like the example shown in Fig. 5. The thermodynam- trate that the classification scheme is generic regardless
ics of category III transitions looks similar to that in 3D of these variations.
free space. In both cases, hR i decreases upon cooling
g
and C peaks at nearly the same temperature, except
V
that a higher peak results for systems of category III. B. Classification of categories using relative surface
SinceadsorptionandH-coreformationnowoccuralmost attraction strengths
together at nearby temperatures, more conformational
degrees of freedom are introduced by the surface inter- Our classification scheme effectively generalized the
actions and this higher entropy gain results in a larger folding behavior into certain transition hierarchies. We
CV. further observe that the dominating factor determining
Details of transitions are again provided by the transition categoryis closely associatedwith the rel-
dhn i/dT, dhn i/dT and dhn i/dT. From ative surface attraction, specifically, the ratio between
HH SH SP
the positions of peaks in dhn i/dT and dhn i/dT, ε +ε and ε .
SH SP SH SP HH
9
TABLEII.CharacteristicthermodynamicsforCategoriesI,II,andIII:Specificheat,CV/N,andaverageradiusofgyrationper
monomer, hRgi/N, as a function of the effectivetemperature kBT/εHH. Errors smaller than thedata points are not shown.
Category I Category II Category III
48mer, εHH =1,εSH =2,εSP =2 48mer, εHH =2,εSH =1,εSP =1 36mer, εHH =3,εSH =1,εSP =1
0.8 0.12 1.2 0.08
1.2 0.08
<R> / N <R> / N
<R> / N g g
g 0.8 0.06 0.8 0.06
0.4 0.08
0.4 Cv / N 0.04 0.4
Cv / N 0.04
Cv / N
0.0 0.04 0.0 0.02 0.0
0.0 2.0 4.0 6.0 0.0 0.4 0.8 1.2 0.0 0.5 1.0
k T / ε k T / ε k T / ε
B HH B HH B HH
67mer, εHH =1,εSH =2,εSP =0 67mer, εHH =1,εSH =0,εSP =1 48mer, εHH =2,εSH =1,εSP =0
0.8 0.12 1.2 0.08
1.2 Cv / N 0.08
Cv / N Cv / N
0.8 0.06
0.8 0.06
0.4 0.08
<R> / N 0.4 0.04 0.4
g <R> / N 0.04
<R> / N g
g
0.0 0.04 0.0 0.02 0.0
0.0 1.0 2.0 3.0 0.0 0.4 0.8 1.2 0.0 0.5 1.0
k T / ε k T / ε k T / ε
B HH B HH B HH
103mer, εHH =1,εSH =1,εSP =1 103mer, εHH =1,εSH =1,εSP =0 103mer, εHH =2,εSH =1,εSP =0
0.8 0.12 1.2 0.08
1.2 0.08
<R> / N
g
<R> / N 0.8 0.06
g 0.8 Cv / N <R> / N 0.06
0.4 0.08 g
Cv / N 0.4 0.04 0.4
Cv / N 0.04
0.0 0.04 0.0 0.02 0.0
0.0 1.0 2.0 3.0 0.0 0.4 0.8 1.2 0.0 0.5 1.0
k T / ε k T / ε k T / ε
B HH B HH B HH
• shows two peaks, one for adsorp- • shows two peaks, one for adsorp- • shows only one peak for a com-
tion at high T,oneforH-corefor- tion at high T,oneforH-corefor- bination of adsorption, flattening
mation at low T mation at low T and H-core formation
CV
• might show a bump for flattening • might have a shoulder or a spike • peak might have one or more
between thetwo peaks forflatteningbelowtheH-corefor- shoulders for processes taking
mation T place at nearby T
• shows a global maximum between • decreases upon cooling, might • keepsdecreasing upon cooling
the two CV peaks, signifying the form a local maximum between • might rise back up due to flatten-
hRgi end of flattening and the start of thetwo CV peaks ing after H-core formation, form-
H-core formation upon cooling
• might rise back up due to flatten- ing a minimum at low T
ing after H-core formation, form-
ing a minimum at low T
10
tween these quantities and the transition categories. We
TABLEIII.Distributionoftransitioncategorieswithrespect
thusbelieve thatε +ε ismoresuitableforthe clas-
totherelativesurfaceattractions. Theabbreviationsreferto SH SP
sification scheme.
thesurfacetypesintroducedinTableI.Thenumbersarethe
We stress that unlike other existing work, this clas-
short forms of the benchmark sequences (e.g. 36 stands for
the 36mer, 48.1 and 48.9 correspond to Seq. 1 and 9 among sification is an inference based on multiple HP se-
theten 48mers, respectively). quences of various chain lengths and attributes. In
Table III, we have also included results from three
εSH other sequences, which were used as a “testing set”
+ Category I Category II Category III
for the adequacy of the classification scheme: a 36mer
εSP
A1: 36, 48.1, (P3H2P2H2P5H7P2H2P4H2P2HP2);another48mer(HP-
48.9, 64, H2P2H4PH3P2H2P2HPH3PHPH2P2H2P3HP8H2); and
67, 103 a 64mer (H12PHPHP2H2P2H2P2HP2H2P2H2P2HP2H2-
A2: 36, 48.9, P2H2P2HPHPH12). The 36mer and the 64mer were
>εHH 64, 67 originally proposed to test a 2D genetic algorithm [26],
H2: 48.1, 48.9, whereas the 48mer is Seq. 1 of the ten testing sequences
64, 67 in [57]. All results from these additional sequences fall
P2: 48.1, 48.9, into the diagonal compartments in Table III, reinforcing
64, 67
that our classification scheme is applicable to other se-
A1/2: 103
quences interacting with an adsorbingsubstrate without
P1: 103
lossofgenerality. This is thus a breakthroughin ourun-
A1/2: 36, 48.1,
derstanding of adsorption properties of lattice proteins:
48.9, 64,
Instead of sequence-dependent individual behavior, the
67
=εHH H1: 48.1, 48.9, thermodynamics of HP proteins do follow common pat-
64, 67, ternsinstructuraltransitionswhentheyinteractwithan
103 adsorbing substrate.
P1: 48.1, 48.9, Even if ε and ε have the same magnitude, they
SH SP
64, 67 are generally not expected to contribute equivalently to
H1/2: 67 H1/2: 48.1, 48.9, thetotalenergybecauseofthecompetitionbetweenε
SH
64, 103
<εHH P1/2: 67, 103 P1/2: 48.1, 48.9 andεHH fortheHmonomers. However,thecontraryap-
A1/3a: 36, 64 pears to be the case in ourstudy because of the entropic
effects, in which the adsorbing P monomers also hinder
a εSH =εSP =1,εHH =3 the formation of H-H contacts in an indirect manner by
dragging the H monomers to the surface. Consider the
surface P2 case for instance: As it is energetically more
favorable to form surface-P contacts than internal H-H
TableIIIshowsthedistributionoftransitioncategories
contacts,thechaintendstositonthesurfaceratherthan
against the relative surface attractions for systems with
to form a compact globule. The H monomers are then
variouschainlengthsandsurfacetypes. Ideally,aperfect
restricted to sit also on the surface instead of forming
correspondencebetweenthetransitioncategoriesandthe
H-H contacts, causing both dhn i/dT and dhn i/dT
relativesurfaceattractionsisimpliedifonlythediagonal SH SP
to peak at the adsorption temperature.
compartmentswerefilled inTable III.In reality,as ther-
modynamic subtleties vary from sequence to sequence,
some off-diagonal compartments are also occupied (e.g.,
VI. CONCLUSION
some systems with ε +ε < ε show category II
SH SP HH
behavior). A few systems also reveal “category duality”
where the thermodynamics bears properties from both Inthis work,proteinadsorptionhasbeenstudiedwith
categories(e.g. the103merinteractingwithsurfacesA1/2 a coarse-grained lattice model, the HP model, interact-
or P1). Nonetheless, the generality of our classification ing with a surface which either attracts all monomers,
scheme is clearly apparentand allows us to infer the fol- only hydrophobic monomers or only polar monomers.
lowing basic rules: transition category I occurs for sur- We have employed Wang-Landau sampling with two ef-
faceswhicharestronglyattractive(ε +ε >≈ε ); fective Monte Carlo trial moves, pull moves and bond-
SH SP HH
categoryIIoccurswhenthe hydrophobicinternalattrac- rebridging moves, to obtain the energy density of states
tionisapproximatelycomparabletothesurfacestrengths and subsequent thermodynamics of structural quantities
(ε +ε ≈ ε ); category III can only occur when for a 48mer, 67mer and 103mer. Ground state energies
SH SP HH
surface strengths are relatively weak compared to ε are also reported. Based on the folding and adsorption
HH
(ε +ε < ε ). Besides, we have also investigated behavior revealed by a careful, comprehensive analysis
SH SP HH
theinfluence oftheproportionofmonomersbyconsider- of the specific heat, radius of gyrationand derivatives of
ing the ratio between (n /N)ε +(n /N)ε and the numbers of surface contacts, we have been able to
SH SH SP SP
ε . However, we did not find an obvious relation be- identify four main types of transition hierarchies (three
HH