Macromolecular
Theory and Simulations

Full Paper

Spherically Symmetric Solvent is Sufficient
to Explain the LCST Mechanism in Polymer
Solutions
Swaminath Bharadwaj, Palakurissi B. Sunil Kumar,* Shigeyuki Komura,Abhijit P. Deshpande
The mechanism of the lower critical solution temperature (LCST) in thermoresponsive
polymer solutions has been studied by means of a coarse-grained single polymer chain simulation and a theoretical approach. The simulation model includes solvent explicitly and thus
accounts for solvent interactions and entropy directly. The theoretical model consists of asingle chain polymer in an implicit solvent where the effect of solvent is included throughthe intrapolymer solvophobic potential proposed by Kolomeisky and Widom. The results ofthis study indicate that the LCST behavior is determinedby the competition between the mean energy difference
between the bulk and bound solvent, and the entropy loss
due to the bound solvent. At low temperatures, solvent molecules are bound to the polymer and the solvophobicity ofthe polymer is screened, resulting in a coiled state. At hightemperatures the entropy loss due to bound solvent offsets
the energy gain due to binding which causes the solvent molecules to unbind, leading to the collapse of the polymer chain
to a globular state. Furthermore, the coarse-grained natureof these models indicates that mean interaction energies are
sufficient to explain LCST in comparison to specific solvent
structural arrangements.

S. Bharadwaj, Prof. A. P. Deshpande
Department of Chemical Engineering
Indian Institute of Technology Madras
Chennai 600036, India
Prof. P. B. S. Kumar
Department of Physics
Indian Institute of Technology Madras
Chennai 600036, India
E-mail: sunil@iitm.ac.in
Prof. S. Komura
Department of Chemistry
Graduate School of Science and Engineering
Tokyo Metropolitan University
Tokyo 192-0397, Japan
Macromol. Theory Simul. 2016, DOI: 10.1002/mats.201600073© 2016 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim

1. Introduction
Thermoresponsive polymers are an important class ofmaterials which exhibit temperature-dependent structural changes and find application in drug delivery,[1,2] surface modification,[3] and self-assembled structures.[4] Ourinterest lies in the family of thermoresponsive polymers
that exhibit a lower critical solution temperature (LCST)
in aqueous solutions. Along with the LCST, these polymers
also exhibit a coil-to-globule transition at the single chain
level. A well-known example of such thermoresponsive
polymers is poly(N-isopropylacrylamide) (PNiPAM) which
exhibits an LCST in water and an upper critical solution

wileyonlinelibrary.com

DOI: 10.1002/mats.201600073

Early View Publication; these are NOT the final page numbers, use DOI for citation !!

1

Macromolecular
Theory and Simulations

S. Bharadwaj et al.

www.mts-journal.detemperature (UCST) in solvents such as methanol, ethanol,
dimethyl sulfoxide, acetone.[5,6] There have been several
experimental[7–9] and simulation[10–15] studies on themechanism of the LCST behavior of PNiPAM in water. The
tunability of the LCST for different applications has also
been explored by studying its variation with additives suchas salt,[16,17] surfactant,[18–20] cosolvents,[21–26] and by thechange in macromolecular architecture such as branching
and tacticity.[27,28]
The origin of LCST in thermoresponsive polymer
solutions is an important question in the field of polymer
science. To understand this phase transition, there
have been several attempts, ranging from the meanfield theory to atomistic molecular dynamic simulations,[10–13,17] focusing on the LCST of PNiPAM-water
system. A mean-field model was proposed by Okada andTanka[29] who hypothesized that preferential interaction among bound water molecules (cooperative hydration) controls the transition. The simulation results of
Deshmukh et al.[13] indicate that the stability of boundwater structure is the driving force for the transition.
Though insightful, one should keep in mind that theresults of these simulations are obtained for a particularpolymer–solvent system, and these approaches requirevery specific interactions, extensive chemical details, orforcefield parameters. Schild et al.[30] have studied the
LCST behavior of PNiPAM in water-alcohol mixtures using
a combination of experiments and the three-component
Flory-Huggins model. In their model, the interaction
parameters have been partially taken from experimental
data, and it does not give us a clear idea about the genericmechanisms that lead to this phenomena.
In the experimental studies, a rich diversity of material
systems based on PNiPAM have been explored to obtainphysical insights related to the effect of substituents,
copolymers, solvents, and additives. From the view pointof the mechanism, Ono and Shikata[9] have calculated thenumber of water molecules per monomer using high-frequency dielectric relaxation measurements. Their results
showed that the LCST is driven by the complete dehydration of the PNiPAM chains, showing the importance ofbound water near the polymer. Bischofberger et al.[6,31]
have performed turbidity and dynamic light scattering on
the ternary system of PNiPAM, water, and alcohol. Theirresults indicate that the thermodynamic description of
the solvent is more important than the specific description of local solvent structure.
Given the multiplicity of systems that can exhibit LCST, it
is highly pertinent to come up with a model that can exhibit
LCST based on broad physical principles. An approach
aimed at identifying the minimal model that exhibits
LCST will help to understand the relative importance ofdifferent contributions. Generic polymer models[32–36]
are suitable candidates for this kind of approach. The

2coarse-grained nature of these models allows us to qualitatively study the importance of the competition between
entropy and internal energy without invoking to a specific polymer or solvent. In this paper, we develop genericpolymer models with spherically symmetric solvent andmonomeric beads for simulation and theoretical studies
of a coil-to-globule transition. Our results indicate that
the LCST depends on the competition between the meaninteraction energy difference between the bulk andbound solvent, and entropy of bound solvent. We showthat a coarse-grained representation of the solvent is sufficient to exhibit an LCST behavior. This indicates that themean interaction energy difference between the boundsolvent and bulk solvent is more important in comparison
to the structural arrangement of the bound solvent. An
important point to note is that this work is aimed at a
generic understanding of the LCST behavior in thermoresponsive polymers, which does not refer to any particularpolymer.
The rest of the paper is organized as follows: in
Section 2, we propose a polymer-solvent model formolecular dynamic simulation studies and introduce thesolvophobic potential[37] used in the theoretical approach.
Section 3 presents the simulation results and numericalcalculations of the theoretical model. In Section 4, ourfindings will be summarized.

2. ModelsIn this section, we discuss our models used for the simulation and theoretical studies. While simulations arecarried out using a bead-spring model for polymers ina homogenous single component solvent, a phenomenological model is used for the theoretical analysis. Below we
describe these models in detail.
2.1. Generic Polymer Model with Explicit Solvent
For the simulation studies, we model the polymer as alinear chain consisting of alternating solvophobic andamphiphilic beads (N total beads, N/2 solvophobic, and N/2
amphiphilic beads). The motivation behind the presenceof two different kinds of beads is to capture the behaviorof the acrylamide family of thermoresponsive polymers
in a generic manner (see Figure 1). The methylene units
along the backbone are analogous to hydrophobic beads.
The substituted methylene units, most generally will haveboth hydrophilic and hydrophobic groups, and therefore
analogous to amphiphilic beads. We emphasize that theintention is to use a generic model, without relating to any
specific polymer system.
The amphiphilic bead has attractive interactions
with both the solvent and the solvophobic beads. The

Macromol. Theory Simul. 2016, DOI: 10.1002/mats.201600073© 2016 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim

www.MaterialsViews.com

Early View Publication; these are NOT the final page numbers, use DOI for citation !!

Macromolecular
Theory and Simulations

Spherically Symmetric Solvent is Sufficient to Explain the LCST Mechanism in Polymer Solutions

www.mts-journal.de

Table 1. Interaction parameters of the SLJ potential. Amphiphilic,
solvophobic, and solvent are represented by A, H, and S, respectively.

ij

AA

HH

SS

AH

HS

AS

−

1

1

1

1

1

1.4, 1.7, 1.8, 2.0

2.5

21/6

2.5

ε ij
rc,ij

Figure 1. Ad hoc mapping of the acrylamide family of thermoresponsive polymers. R1 and R2 can be any arbitrary groups. The partof the monomer within the red box may have both hydrophilic
and hydrophobic groups due to which it is modeled as an amphiphilic bead.

interaction between the solvophobic bead and the solvent is purely repulsive. The solvent can be introduced
either by including it explicitly or by incorporating its
effects implicitly within the interaction potentials. Sincecoarse-grained potentials are obtained by integrating out
the internal degrees of freedom of the unit, they are temperature dependent in general. When the scale of coarse
graining is small, such a dependence is weak and canbe neglected. However, when the solvent is implicit, theinteraction potentials have a stronger dependence on thetemperature, and the nature of this dependence has to be
assumed a priori. To avoid a specifically assumed temperature dependence of the interaction potentials, we explicitly incorporate the solvent. The potential energy for thesystem is given by Equation (1)
12

6

12

6

N −1
Nt
 σ σ   σ 
σ  
E = ∑ kb (bi − bi 0 )2 + ∑ ∑ 4 ε ij   −   −   +   
r
r
r






 rc,ij  

ij
ij
c,
ij
i= 1
i= 1 j > i




(1)

where N is the number of beads in the polymer chain,
as mentioned before. kb the force constant for thebonded interaction, bi the bond length between neighboring beads, Nt the total number of beads in the system(polymer + solvent), bi0 the equilibrium bond length,
and rij the distance between two nonbonded beads. The
second term is the Shifted Lennard Jones (SLJ) potential
with rc,ij being the cutoff distance at which the potential
is truncated and shifted to zero. The above form of SLJ
potential ensures that all the beads are spherically symmetric and have size σ. All the interaction parameters are
kept independent of the temperature. We define dimen−
sionless quantities as rij = rij /σ , ε ij = εij /ε ss, kb = σ 2kb /ε ss,3
bi 0 = bi 0 /σ , T = kBT /ε ss, P = σ P/ε ss, and t = t ε ss /(mσ 2 ) ,where εss is the potential energy of interaction between
two solvent beads. We fix the values to bi0 = 1 and kb = 200
for all the simulations. The values of the other interaction
parameters are listed in Table 1, where εAS, εAH, and εAA
are the interaction energies of the amphiphilic bead with

www.MaterialsViews.com

2.5

2.5

2.5

the solvent bead, the solvophobic bead, and the amphiphilic bead, respectively.
Molecular dynamic simulations were performed in
an NPT ensemble using the Nose-Hoover thermostat fordifferent temperatures at a constant pressure P = 0.002.
The trajectories were generated using the VelocityVerlet algorithm with a time-step ∆t = 0.004 . The ratioof the number of polymer beads to the solvent beads
was maintained at 0.04 for all simulations. Simulations of N = 200 chain were performed at four different
interactions; ε− AS = 1.4,1.7,1.8,2.0. For each of these values,
the temperature was varied from T = 0.5 to 0.7 with aninterval of 0.05. Simulations were also performed for an
N = 400 chain for ε− AS = 1.7 at the temperatures ranging
from T = 0.5 to 0.8 with an interval of 0.05. The 200 and
400 bead systems were equilibrated for 1 × 108 steps, andthe data were sampled after every 4 × 106 and 2 × 107
steps, respectively. Four different initial configurations
were used for averaging. All simulations were performed
using open source molecular dynamics code LAMMPS.[38]The simulation data were used for the calculation ofdifferent structural quantities. We calculated the radiusof gyration, Rg, of the polymer to monitor the swelling of
the polymer chain. The average Rg and standard deviation
were calculated from the distribution obtained by sampling 1000 simulation replicas. The standard deviation
was plotted as the error bar. We define a dimensionlessradius of gyration Rg = Rg /σ given by Equation (2)
Rg =

1 N
∑(ri − rcm )2
N i= 1(2)where N = 200 and 400, rcm and ri are the dimensionlesscoordinates of the center of mass of the polymer chain
and the ith bead, respectively.
The number of solvent beads, Ns, in the first solvation
shell of the polymer was calculated to determine thebound solvent content. The solvent beads which were
within a distance of r = r / σ = 1.5 from any of the polymeric beads were regarded to be part of the first solvation
shell.
The effective interaction between the polymer beads
was calculated using the potential of mean force UAH,
which was calculated from the radial distribution function of the amphiphilic and solvophobic bead pairs using
Equation (3)[39]

Macromol. Theory Simul. 2016, DOI: 10.1002/mats.201600073© 2016 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim

Early View Publication; these are NOT the final page numbers, use DOI for citation !!

3

Macromolecular
Theory and Simulations

S. Bharadwaj et al.

www.mts-journal.de

U AH (r )= − ln gAH (r )T

(3)where U AH = U AH /ε ss is the dimensionless potential ofmean force.
2.2. Solvophobic Potential by Kolomeisky and Widom
Our simulations indicate that the LCST is dependent
on the entropy loss of the bound solvent, and the meaninteraction energy difference between the bulk and thebound solvent (see Section 3.1). To obtain further insights
related to the nature of the transition, scaling behavior,and the effect of chain flexibility, a theoretical approachis adopted where we consider the hydrophobic potential proposed by Kolomeisky and Widom (KW).[37] The
KW model is spherically symmetric in nature and does
not relate to a specific solvent or solute, which is in correspondence with the modeling framework employedin our simulations. The KW model belongs to a class of
implicit solvent models[40,41] which incorporate different
interaction energies depending on the proximity of thesolvent to a solute molecule. These models have beenused for studying the solubility of small solutes in water.
The KW model is one of the simplest among these models
as it has only two solvent interaction energies; one for thebound state and another for the bulk state. These interaction energies are analogous to the monomer-solvent
interaction energy (εAS), and solvent–solvent interactionenergy (εSS) in our simulation model, respectively. Henceit can be seen that the model contains those contributions which have been emphasized in our simulation
results.
In the KW model, solvent molecules form a 1D latticewith a nearest neighbor interaction, and each solvent
molecule can exist in q different states denoted by 1, 2, ⋅⋅⋅qas shown in Figure 2. The interaction energy between theneighboring solvent molecules is w when both of them
are in the state “1”, and u otherwise with u > w. Herethe former and the latter cases are termed as the boundstate (BS) and the unbound state (US), respectively. A solvent molecule in the BS state can exist only in one state,
whereas that in the US state in q − 1 states. Hence theentropy of a solvent molecule in the BS state vanishes,
whereas that in the US state is given by kBln (q − 1). In
other words, the BS state is energetically favorable (w < u),while the US state is entropically favorable. The energetic(ΔU) and entropic (ΔS) differences between the BS and
US states are w − u and −kBln (q − 1), respectively. The
competition between the BS and the US states can beconveniently described by a dimensionless parameter x
defined as

4

Figure 2. Schematic representation for the theoretical model.
Black and blue beads represent monomer and solvent molecules,
respectively. a) Single polymer chain in explicit solvent, b) single
polymer chain in implicit solvent where the effect of solvent isincorporated in the monomer–monomer interaction potential
φ(r), and c) φ(r) modeled by the solvophobic potential given by
KW.[37] The potential assumes a 1D solvent lattice where eachsolvent molecule has q states. Neighboring solvent molecules
exist in a bounded state (BS) when both are in the state “1” andin unbounded state (US) state otherwise. Solute can occupy interstitial sites between bounded (BS) solvent molecules.

x = e( ∆U −T ∆S )/kBT =

q−1
c

(4)

where
(5)

c = e( u−w )/kBT

Since the number of states q is constant in the KW
model, x is a monotonically increasing function of thetemperature. Solute molecules are allowed to occupy only
the interstitial sites between the solvent molecules in the
BS state. Based on these assumptions, KW obtained thesolvent mediated attraction potential φ(r) between two
solute molecules (implicit solvent) for r > σ[37]  1 + Q   1 − S  ( r −σ )/σ 
φ (r ) = −kBT ln  1 +   1 − Q   1 + S 


(6)

where
4x  1  
S = 1 −1− x )2  c  
(1
+

Q=

sgn( x -1)
 1 + 4 x / ( x − 1)2 c 

1/2

(7)

(8)

1/2

and sgn( z ) is the sign function. Figure 3 shows thevariation of φ/kBT as a function of r = r / σ for two different temperatures. The range of the solvent mediated interaction becomes shorter when the temperature is increased. From the inset of Figure 3, it can beseen that the attraction becomes stronger for higher

Macromol. Theory Simul. 2016, DOI: 10.1002/mats.201600073© 2016 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim

www.MaterialsViews.com

Early View Publication; these are NOT the final page numbers, use DOI for citation !!

Macromolecular
Theory and Simulations

Spherically Symmetric Solvent is Sufficient to Explain the LCST Mechanism in Polymer Solutions

www.mts-journal.de

3. Results and Discussions

0

3.1. Simulation

φ/kBT

-3
−11

-6

−12
−13
1.0

-9

1.5

2.0

T̃ = 0.080
T̃ = 0.087

-12
1

8

15

22

29

r̃
Figure 3. Variation of φ/kBT as function of r = r/σ for differentT (= T /(u − w /kB )) when q = 5 × 105. The inset shows the variation
at low r values.

temperature when r is small. Such a behavior indicatesthat the monomers tend to aggregate as the temperature is increased. The large value of q used by KW was
justified in order to match the temperature dependenceof the solubility of nonpolar solutes in water.[37] From
Equation (6) it can be seen that the monomer–monomer
potential is temperature dependent due to the implicit
nature of solvent. An important point to note is that
this temperature dependence is not ad hoc but a consequence of the underlying solvent model summarized
in Figure 2.

6.5 (a)AS = 1.4

AS = 1.7

AS = 1.8

AS = 2.0

Simulations on N = 200 chain were performed to explore
the structural change of the polymer chain with the temperature for different εAS values. This sweep across εASvalues was performed to identify interesting regimes.
Figure 4a shows the variation of Rg with T for different
ε− AS . The error bars are larger for the expanded state (higher
εAS values) as the number of conformations are higher due
to the absence of angular and dihedral interactions. However, the measured Rg is distinctly larger for larger ε− ASvalues, implying swelling as a result of stronger association between the amphiphilic bead and the solvent. As afunction of temperature, we observe that there are threedifferent behaviors according to the value of ε− AS; (i) whenε− AS = 1.4, the polymer chain remains in a collapsed stateand its Rg shows a slight increasing trend with T and theerror bars are very small, (ii) when ε− AS = 2.0, the polymer
chain is in the large Rg state with large error bars at all thetemperatures, and (iii) when ε− AS = 1.7 and 1.8, Rg shows aslight decrease with increase in temperature. Though thisdecrease of Rg with temperature for εAS = 1.7, 1.8 is very
interesting, it is inconclusive and requires further probing
due to the large value of error bars.
To further explore these interesting trends, simulations
were performed on a larger system (N = 400 chain). We
choose εAS = 1.7 for N = 400 chain as it has the most significant decrease in average Rg with temperature among
the studied εAS values. From Figure 4b, we see that thesize change with the temperature for the N = 400 chain
is more prominent compared to that of the N = 200 chain,
which is due to difference in scaling of Rg with N for the

8.5
7.5

Rg

Rg

5.5
4.5

6.5
5.5

3.5
2.5
0.50

AS = 1.7

(b)

4.5

0.55

0.60

0.65

0.70

0.5

0.6

T

T
−

Figure 4. Variation of Rg of polymer chain with T for a) 200 bead polymer at different values of ε AS , the error bars for εAS = 1.4 are the size−of the symbol. b) N = 400 chain at ε AS = 1.7.

www.MaterialsViews.com

Macromol. Theory Simul. 2016, DOI: 10.1002/mats.201600073© 2016 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim

Early View Publication; these are NOT the final page numbers, use DOI for citation !!

5

Macromolecular
Theory and Simulations

S. Bharadwaj et al.

www.mts-journal.deFigure 5. Representative snapshots of equilibrated N = 400 chain at εAS = 1.7 for different temperatures. Green and red beads represent
solvophobic and amphiphilic units, respectively. a) T = 0.55, b) T = 0.6 , c) T = 0.7 and d) T = 0.8. Solvent beads are not included for clarity.

coiled and globular states.[42] In Figure 4b, Rg decreases
by 23% as the temperature is increased from T = 0.55 to
0.75. Such trends have been observed in atomistic simulations of PNiPAM-water system.[12,17,43] Unlike the firstorder LCST transition behavior observed in experimental
studies of thermoresponsive polymer solutions,[8,44,45] thetransition observed in Figure 4 b shows a gradual change
in Rg with T . It is known that the coil-to-globule transition is a first-order transition only in the case of rigid andsemirigid chains.[46–48] Hence the continuous transitionin our simulation is not surprising because the polymer
chain used in the simulations is fully flexible due to theabsence of the angular and dihedral interactions. This isfurther supported by the results of our theoretical model,
where we observe that the coil-to-globule transition deviates from a first order behavior with increase in chain
flexibility (see Section 3.2). Based on the above observa−
tions, we consider that the behavior observed for ε AS = 1.7
is akin to the LCST phenomenon in thermoresponsive
polymers. Figure 5 shows representative snapshots of
N = 400 chain at εAS = 1.7 for different temperatures, we cansee that the polymer chain is in the coil-like state below
T = 0.6, while it is the globule-like state above T = 0.7 .
Hereafter we will be referring to ε− AS = 1.4 , ε− AS = 1.7 and
1.8, and ε− AS = 2.0 as low, intermediate, and high values,

400

(a)

AS = 1.4

AS = 1.7

AS = 1.8

AS = 2.0

respectively. To further understand the behavior in thesethree states, we examine the variation of the number ofbound solvent beads, Ns, with the temperature. In Figure 6a,
the number of bound solvent is found to be almost independent of T for low and high values of ε− AS, whereas it
decreases by increasing the temperature for intermediate
values. In Figure 6b, we plot Ns for the N = 400 chain as afunction of T , when εAS = 1.7. The decrease of Ns with thetemperature is more prominent as compared with that of
the N = 200 chain case. Moreover, Ns markedly decreases
around T = 0.65 which coincides with the temperaturearound which Rg also decreases (see Figure 4b). It should
be stressed that the change in Rg and Ns with the temperature is observed in our model with minimal interactions asin Equation (1). It is important to note that even when theinteraction parameters are independent of the temperature, the temperature dependence of Rg is induced by thebound solvent number Ns for a range of values of ε− AS.
The above trends in the three different states can berationalized by considering the different energy andentropic contributions in the model system. The attraction between the amphiphilic monomer and the solvent
is energetically favorable since it is stronger than thesolvent–solvent interaction (εAS > εSS), leading to a coil-like
polymer conformation with many bound solvent beads.

550

350

500

Ns

Ns

AS = 1.7

(b)

300

450

250
200
0.50

0.55

0.60

0.65

0.70

400

0.5

0.6

T

0.

.8

T

Figure 6. Variation of Ns with T at different values for a) N = 200 chain at different values of ε AS, b) N = 400 chain at ε AS = 1.7.
−

6

Macromol. Theory Simul. 2016, DOI: 10.1002/mats.201600073© 2016 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim−

www.MaterialsViews.com

Early View Publication; these are NOT the final page numbers, use DOI for citation !!

Macromolecular
Theory and Simulations

Spherically Symmetric Solvent is Sufficient to Explain the LCST Mechanism in Polymer Solutions

www.mts-journal.de

-3.1

-2.75

-3.00

-3.25

-3.50

1.00 1.05 1.10 1.15 1.20 1.25

T = 0.5
T = 0.6
T = 0.7(b)

U AH/kBTU AH/kBT(a)

T = 0.5
T = 0.6
T = 0.7

-3.3

-3.5

-3.7
1.00

1.05

1.10

1.15

1.20

r

r

Figure 7. Variation of UAH / kBT with r at different temperatures for a) N = 400 chain at ε− AS = 1.7, b) N = 200 chain at ε− AS = 1.4.

Whereas from the viewpoint of entropy, the solvent prefers to be in the bulk state rather than the bound state.
Another important contribution is the solvation of solvophobic beads, which is unfavorable as the interaction
between the solvophobic bead and the solvent is repulsive. The interplay of these contributions determines thedifferent states of the polymer chain.
Concerning intermediate ε− AS values, the large values of
Rg in the coil-like state (see Figure 4b) and large Ns at low
temperatures (see Figure 6b) indicate that the attraction
between the amphiphilic monomer and the solvent dominates the entropy loss of the solvent. Hence it is favorable
for the solvent to solvate the solvophobic beads despite
the repulsive interaction between them. With increase in
the temperature, the gain in the entropy due to unbindingof the solvent leads to a reduction in the bound solvent
content, Ns. In the transition region (0.6 < T < 0.7), there
is marked decrease of the bound solvent because theentropy gain of the solvent dominates over the amphiphilic-solvent attraction. This change in the dominant
contribution makes the solvation of solvophobic beads
unfavorable, leading to an attraction between the polymeric beads. Such an attraction becomes stronger withthe temperature, and drives the transition from the coillike state to the globule-like state. In Figure 7a, we plot
the potential of mean force, U AH, between the amphiphilic and the solvophobic beads as a function of the
distance r . Here we see that the attraction between thepolymeric beads increases with the temperature. At low
temperature, the attraction between the polymeric beads
is screened due to the presence of solvent beads. At hightemperatures, the bound solvent beads unbind due toentropy gain which reduces the screening effect, leading
to the collapse of the polymer. As shown in Figure 7b, theincrease in the attraction between polymeric beads withtemperature is not observed for low ε− AS values.
For low ε− AS values, both Rg and Ns are small. Thismeans that at these interaction strengths the entropic

www.MaterialsViews.com

gain of the free solvent beads dominates over theinteraction between the amphiphilic monomer and thesolvent for all the temperatures. Hence the solvation ofsolvophobic beads is unfavorable in the entire temperature range. On the other hand, for high ε− AS values, both
Rg and Ns are large. This indicates that the attractive
interaction is stronger than the entropy loss which leads
to the binding of the solvent to the polymer. Hence, thesolvation of the solvophobic beads is favored at all temperatures despite the repulsive interaction between thesolvophobic monomer and the solvent.
We find from our simulations that the transition fromthe coil-like state to the globule-like state depends on
how much the solvophobic beads can be solvated. The
solvation depends on the interplay between the entropy
loss of the bound solvent, and the energetic difference
between the bound and the bulk solvent. The boundsolvent is energetically favored as the interaction between
the amphiphilic monomer and the solvent is stronger
than the solvent–solvent interaction. Therefore two kindsof beads having opposite interactions with the solvent arenecessary to exhibit the LCST as long as the interaction
parameters are independent of the temperature as in ourcase. Additionally in our simulation, we did not incorporate any specific chemical or structural details pertaining
to the solvent or the polymer chain. Therefore we havedemonstrated that the LCST behavior can be exhibitedby a coarse-grained description of the polymer and thesolvent even when the interaction parameters are independent of the temperature.
3.2. Theoretical Description of the Coil-to-Globule
Transition
To study the coil-to-globule transition in the theoreticalframework, we adopt the phenomenological free energy
expression for a polymer chain in an implicit solvent,which has been given by Grosberg and Kuznetsov.[46] We

Macromol. Theory Simul. 2016, DOI: 10.1002/mats.201600073© 2016 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim

Early View Publication; these are NOT the final page numbers, use DOI for citation !!

7

Macromolecular
Theory and Simulations

S. Bharadwaj et al.

www.mts-journal.dechoose the KW potential (see Section 2.2) to model theeffective monomer–monomer interaction. The free energy
expression for the system is as follows

1.0

F
1
N
1
=α2 + 2 + 3 3 B + 6 6 C
kBT
α ασ
α σ

(9)

(

)

(10)

where Φ(r) is the monomer–monomer interaction potential. Here we assume that it has the following form
 ∞ : r < σ
Φ(r ) = φ (r ) : r > σ


(11)

where Φ(r) for r <σ corresponds to the (hard-core) excluded
volume interaction, and φ(r) for r > σ (see Equation (6)) is
the KW solvophobic potential. Substituting Equations (6)and (11) into Equation (10), we obtain the second virial
coefficient B as
 1+Q  1
2πσ 
2
2 
B=
−
1 + 3
2 +
 1 − Q   ln L ( ln L ) ( ln L )3  
3 
3

-0.5
-1.0
0.060

q = 5 x 105
q = 1 x 106
q = 2 x 106
0.064

0.068

0.072

0.076

T̃
Figure 8. Variation of B with T for different q values.

α 5 − α − π N B −

C
=0
α3

(13)

where C = 3C /σ 6 is the nondimensional third virial coefficient which is related to the rigidity of the polymer
chain and has contributions from the three-body interactions such as angular interactions. We numerically solve
the above equation to obtain α for different values of B

1.6

(12)

where L = (1 − S)/(1 + S).
Hereafter we use the dimensionless quantities suchas the temperature T = kBT /(u − w ), the distance r = r/σ ,the solvent mediated interaction potential φ = φ /(u − w ),and the second virial coefficient B = 3 B /2πσ 3 . In Figure 8,
we plot B as a function of T for different q-values. We
observe that B remains almost unity for low temperatures and then decreases rapidly for higher temperatures.
This means that the monomer–monomer interaction isrepulsive for lower temperatures, while it is attractive forhigher temperatures. Moreover, the temperature corresponding to the sharp drop from the positive (repulsive)to the negative (attractive) B-values
decreases when q isincreased for a fixed u − w value. The significance of thisobservation will be discussed later.
In order to find the equilibrium polymer conformation,
we minimize Equation (9) with respect to α, and obtain
the equation

8

0.0

1.2

α

∞

0

0.5

B̃

where α = Rg / ( N σ ) characterizes the extent of theswelling, N the degree of polymerization, σ the diameter ofthe monomeric unit, B and C are the second and the third
virial coefficients, respectively. The first and the second
terms in the above free energy are the entropic contributions with the coiled and the globular states given by
α > 1 and α < 1, respectively. The third and the fourth terms
represent the energy contributions from the two-body andthe three-body interactions, respectively. In general, thesecond virial coefficient B is given by
B = 2π ∫ dr r 2 1 − e − Φ( r )/kBT

1.5

0.8
0.4

C̃ = 1C̃ = 10−1C̃ = 10−2C̃ = 10−40.0
756.081 756.083 756.085 756.087

T̃ × 104
Figure 9. Variation of α with T for different values of C whenq = 5 × 105 and N = 105.

Macromol. Theory Simul. 2016, DOI: 10.1002/mats.201600073© 2016 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim

www.MaterialsViews.com

Early View Publication; these are NOT the final page numbers, use DOI for citation !!

Spherically Symmetric Solvent is Sufficient to Explain the LCST Mechanism in Polymer Solutions

Macromolecular
Theory and Simulations
www.mts-journal.de In Figure 9, we show theat fixed values of q, N and C.
variation of the swelling parameter α with the temperature T for different values of the third virial coefficientC when q = 5 × 105. We see that the polymer chain undergoes a transition from the coiled state to the globularstate with increase in the temperature.
For low temperatures, the polymer chain is in the
coiled state and hence large α. In this case, the first andthe third terms in Equation (13) are dominant because
N is also large. Then Equation (13) reduces toα 5 − π N B ≈ 0

(14)

∼N3/5 corresponding to the scaling of aand we obtain Rg
polymer chain in a good solvent. When the temperature
is high, on the other hand, the polymer is in the globularstate and α becomes small. Then the third and the fourthterms dominate in Equation (13)

π N B +

C
≈0
α3

(15)

This gives the scaling Rg ∼ N1/3 corresponding to apolymer chain in a poor solvent.
In order to understand this coil-to-globule transitionin terms of the dominant interactions, let us consider thefree energy difference of a solvent molecule between the
BS and US states∆F = FBS − FUS = U BS − U US − T ( SBS − SUS )(16)

where UBS (UUS) and SBS (SUS) are the energy and theentropy of the BS (US) state, respectively. Using the model
parameters defined in Section 2.2, the above quantity canbe expressed as∆F = −(u − w ) + kBT ln( q − 1) = kBT ln x

(17)where x is defined before in Equation (4). From
Equation (17) and since u > w, we find ΔF < 0 for x < 1(lower temperatures), leading to the BS state being morefavorable. Given the implicit nature of the solvent in thetheoretical treatment, the large value of α for lower temperatures in Figure 9 is an indication of the large amountof the bound solvent. For x > 1 (higher temperatures), on
the other hand, the US state is more favorable and theamount of the unbound solvent increases. In this case,the attraction between the monomeric units is induced.
This is seen in Figure 8 for the larger negative value of the These solvent induced intersecond virial coefficient B.
actions drive the transition from the coiled state to theglobular state. The phase transition temperature T* determined by the condition ΔF = 0 is given by Equation (18)
T∗ =

u−w
kB ln( q − 1)

www.MaterialsViews.com(18)

As observed in Section 2.2, T* decreases as q is increased
when u −w is fixed.
We further observe in Figure 9 that the nature of thetransition changes from a discontinuous transition to acontinuous one by increase in the third virial coefficient C .
It is known that the third virial coefficient is larger for
flexible chains.[46,47] This shows that the nature of transition deviates from a first-order behavior with increasein the flexibility of the chain. Another important pointis that the above argument does not include any details
regarding the structure of the solvent and/or polymer
which indicates that the structural details of the boundsolvent are not necessary for the thermoresponsive
behavior to manifest.

4. Conclusions
In this paper, we have tried to understand the single chain
coil-to-globule transition of a thermoresponsive polymer
through simulation and theoretical approaches. In the simulations, the model comprises of a single polymer chain in
an explicit solvent with temperature independent interaction parameters. The solvent is explicitly included to avoid
any ad hoc dependence of the interaction potentials on
temperature. To obtain further insights, we have adopted
a theoretical framework with only those interactions thathave been emphasized in our simulation studies. The
theoretical model describes a single chain in an implicit
solvent where the effect of solvent is included into themonomer-monomer interaction potential. For the interaction between the monomers, we have used the solvophobic potential proposed by Kolomeisky and Widom.[37]The temperature dependence of the solvophobic potential
is not ad-hoc in nature and arises due to the underlying
solvent model.
Our simulations indicate that the LCST is dependent
on the competition between the two contributions,
namely, the entropy loss of the bound solvent and themean energy difference between the bound and thebulk solvent. This hypothesis is supported by our theoretical calculations. The former favors the globular stateof the polymer chain, whereas the latter prefers the coiled
state. At low temperatures, solvent molecules bind to thepolymer chain rather than to reside in the bulk because
the bound state is energetically favorable, and the coiled
state is obtained. At high temperatures, on the other
hand, the entropy loss of the bound state is more dominant than its energetic gain and the solvent molecules
tend to leave the polymer. Such a competition of the interactions induces a solvent driven attraction between thepolymeric beads leading to the collapse of the polymer
chain. These findings are supported by the experimental
studies of Bischofberger et al. who showed that the LCST

Macromol. Theory Simul. 2016, DOI: 10.1002/mats.201600073© 2016 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim

Early View Publication; these are NOT the final page numbers, use DOI for citation !!

9Macromolecular
Theory and Simulations

S. Bharadwaj et al.

www.mts-journal.de

is dependent on the mean energy difference between thebound and bulk solvent.[6,31]
Another common feature between the simulation andthe theory is that both of them have spherically symmetric solvent and monomeric beads. This indicatesthat the structural arrangement of the solvent moleculesaround the polymer chain is not a necessary component
to be considered explicitly for the coil-to-globule transition. In other words, a coarse-grained description of
the solvent is sufficient to reproduce the LCST behavior,
this observation is in agreement with the experimental
results of Bischofberger et al.,[6,31] where they show thata coarse-grained representation of the solvent is sufficient to explain the LCST and its variation with different
alcohols.
In our simulations studies, it is shown that the LCST isdependent on the amphiphilic-solvent attraction and solvation of the solvophobic beads. This indicates that in thecase of temperature-independent interaction para­meters,
two different kinds of beads with opposite interaction
with the solvent are required in the monomeric unit
to exhibit LCST. In the theoretical model, the solvent is
implicit and the variation of solvent effects with temperature is included in the monomer-monomer interaction
due to which one kind of bead is sufficient to exhibit LCST.
A mixture of PNiPAM and water exhibits a first-orderphase transition at 32 °C.[7,8] In our theoretical argument,
we showed that the order of the transition changes from afirst-order (discontinuous) to a second-order (continuous)with the increase in the chain flexibility (or the third
virial coefficient). The reason why we observe only thecontinuous coil-to-globule transition in our simulation is
because we did not take into account any angular interactions to deal with the chain flexibility. The effect ofangular and dihedral interactions on the nature of thecoil-to-globule transition needs to be explored further.
The results of our model are applicable to cross-linked
thermoresponsive polymers as the lengths of the segments between successive cross-links are usually very
long. Also, it is known that the degree of cross linking
only affects the extent of swelling and does not change
the LCST.[49] Another interesting feature of our model isthat the variation in behavior with temperature for different ε− AS values (see Figure 4a) is analogous to different
polymers in the family of poly(N,N-alkyl alkyl acrylamide). Low ε− AS values are similar to polymers with bulky
side groups such as poly(N-butyl acrylamide) which arealways insoluble in water. On the other hand, high ε− ASvalues are similar to polymers such as polyacrylamide
which are always soluble in water. This indicates that
our coarse-grained model is able to explain the LCST
behavior for different systems in the poly(N,N-alkyl,alkyl
acrylamide) family of polymers. Furthermore, the model
can also be utilized to study the variation of LCST due to10

factors such as co-solvents and additives by examining
their effect on the two dominant physical interactions
namely, the entropy loss of the bound solvent and themean energy difference between the bound and thebulk solvent. This aspect will be addressed in our futurestudies.
Acknowledgements: The computations were carried out atthe High Performance Computing Facility at IIT Madras. S.K.
acknowledges support from the Grant-in-Aid for Scientific
Research on Innovative Areas “Fluctuation and Structure” (Grant
No. 25103010) from the Ministry of Education, Culture, Sports,
Science, and Technology of Japan, the Grant-in-Aid for Scientific
Research (C) (Grant No. 15K05250) from the Japan Society for the
Promotion of Science (JSPS), and the JSPS Core-to-Core Program“International Research Network for Non-equilibrium Dynamics
of Soft Matter.” S.B. thanks Tokyo Metropolitan University forthe support provided through the co-tutorial program. S.B.
thanks Ryuichi Okamoto, David Andelman, Kenichiro Koga andKazuhiko Seki for valuable discussions.
Received: August 9, 2016; Revised: October
Published online: ; DOI: 10.1002/mats.201600073

19,

2016;

Keywords: coarse-grained models; coil-to-gloule transition;
lower critical solution temperature; stimuli-sensitive polymers;swelling; thermoresponsive polymer solutions[1] T. Hoare, J. Santamaria, G. F. Goya, S. Irusta, D. Lin, S. Lau,
R. Padera, R. Langer, D. S. Kohane, Nano Lett. 2009, 9, 3651.
[2] T. Okano, A. Kikuchi, Y. Sakurai, Y. Takei, N. Ogata, J. Controlled Release 1995, 36, 125.
[3] T. R. Hoare, D. S. Kohane, Polymer 2008, 49, 1993.
[4] S. Sun, P. Wu, W. Zhang, W. Zhang, X. Zhu, Soft Matter 2013,
9, 1807.
[5] R. O. Costa, R. F. Freitas, Polymer 2002, 43, 5879.
[6] I. Bischofberger, D. C. E. Calzolari, V. Trappe, Soft Matter
2014, 10, 8288.
[7] X. Wang, X. Qiu, C. Wu, Macromolecules 1998, 31, 2972.
[8] C. Wu, X. Wang, Phys. Rev. Lett. 1998, 80, 4092.
[9] Y. Ono, T. Shikata, J. Am. Chem. Soc. 2006, 128, 10030.
[10] S. Deshmukh, D. A. Mooney, T. McDermott, S. Kulkarni,
J. M. Don Mac Elroy, Soft Matter 2009, 5, 1514.
[11] S. Deshmukh, D. A. Mooney, J. M. MacElroy, Mol. Simul. 2011,
37, 846.
[12] M. Alaghemandi, E. Spohr, Macromol. Theory Simul. 2012,
21, 106.
[13] S. A. Deshmukh, S. K. R. S. Sankaranarayanan, K. Suthar,
D. C. Mancini, J. Phys. Chem. B 2012, 116, 2651.
[14] A. K. Tucker, M. J. Stevens, Macromolecules 2012, 45, 6697.
[15] E. Chiessi, A. Lonardi, G. Paradossi, J. Phys. Chem. B 2010,
114, 8301.
[16] Y. Zhang, S. Furyk, D. E. Bergbreiter, P. S. Cremer, J. Am. Chem.
Soc. 2005, 127, 14505.
[17] H. Du, R. Wickramasinghe, X. Qian, J. Phys. Chem. B 2010,
114, 16594.
[18] H. G. Schild, D. A. Tirrell, Langmuir 1991, 7, 665.
[19] V. S. Shinde, M. V. Badiger, A. K. Lele, R. A. Mashelkar, Langmuir 2001, 17, 2585.

Macromol. Theory Simul. 2016, DOI: 10.1002/mats.201600073© 2016 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim

www.MaterialsViews.com

Early View Publication; these are NOT the final page numbers, use DOI for citation !!

Spherically Symmetric Solvent is Sufficient to Explain the LCST Mechanism in Polymer Solutions

Macromolecular
Theory and Simulations
www.mts-journal.de[20] Y. Mohan, T. Premkumar, D. Joseph, K. Geckeler, React. Funct.
Polym. 2007, 67, 844.
[21] X. Zhang, Y. Yang, T. Chung, Langmuir 2002, 18, 2538.
[22] C. Scherzinger, A. Schwarz, A. Bardow, K. Leonhard,
W. Richtering, Curr. Opin. Colloid Interface Sci. 2014, 19,84.
[23] D. Mukherji, K. Kremer, Macromolecules 2013, 46, 9158.
[24] D. Mukherji, C. M. Marques, K. Kremer, Nat. Commun. 2014,
5, 4882.
[25] D. Mukherji, C. M. Marques, T. Stuehn, K. Kremer, J. Chem.
Phys. 2015, 142, 114903.
[26] J. Walter, V. Ermatchkov, J. Vrabec, H. Hasse, Fluid Phase
Equilib. 2010, 296, 164.
[27] B. Ray, Y. Okamoto, M. Kamigaito, M. Sawamoto, K.-I. Seno,
S. Kanaoka, S. Aoshima, Polym. J. 2005, 37, 234.
[28] Y. Katsumoto, N. Kubosaki, Macromolecules 2008, 41, 5955.
[29] Y. Okada, F. Tanaka, Macromolecules 2005, 38, 4465.
[30] H. G. Schild, M. Muthukumar, D. A. Tirrell, Macromolecules
1991, 24, 948.
[31] I. Bischofberger, D. C. E. Calzolari, P. De Los Rios, I. Jelezarov,
V. Trappe, Sci. Rep. 2014, 4, 4377.
[32] J. A. Anderson, A. Travesset, Macromolecules 2006, 39, 5143.
[33] S. J. Marrink, H. J. Risselada, S. Yefimov, D. P. Tieleman,A. H. de Vries, J. Phys. Chem. B 2007, 111, 7812.

www.MaterialsViews.com[34] M. Mella, L. Izzo, Polymer 2010, 51, 3582.
[35] J. M. Polson, N. E. Moore, J. Chem. Phys. 2005, 122, 024905.
[36] M. Hatakeyama, R. Faller, Phys. Chem. Chem. Phys. 2007, 9,4662.
[37] A. B. Kolomeisky, B. Widom, Faraday Discuss. 1999, 112, 81.
[38] S. Plimpton, J. Comput. Phys. 1995, 117, 1.
[39] D. Chandler, Introduction to Modern Statistical Mechanics,
Oxford University Press, Oxford 1987.
[40] B. Lee, G. Graziano, J. Am. Chem. Soc. 1996, 118, 5163.
[41] S. Moelbert, P. D. L. Rios, Macromolecules 2003, 36, 5845.
[42] A. Y. Grosberg, A. R. Khokhlov, Giant Molecules: Here, There,
and Everywhere, Academic Press, New York, USA 1997.
[43] S. A. Deshmukh, S. K. R. S. Sankaranarayanan, K. Suthar,
D. C. Mancini, J. Phys. Chem. B 2012, 116, 2651.
[44] X. Wang, C. Wu, Macromolecules 1999, 32, 4299.
[45] Z. Cao, W. Liu, P. Gao, K. Yao, H. Li, G. Wang, Polymer 2005, 46,5268.
[46] A. Y. Grosberg, D. V. Kuznetsov, Macromolecules 1992, 25,
1970.
[47] G. Graziano, Int. J. Biol. Macromol. 2000, 27, 89.
[48] B. M. Baysal, F. E. Karasz, Macromol. Theory Simul. 2003, 12,627.
[49] T. Çaykara, S. Kiper, G. Demirel, J. Appl. Polym. Sci. 2006, 101,
1756.

Macromol. Theory Simul. 2016, DOI: 10.1002/mats.201600073© 2016 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim

Early View Publication; these are NOT the final page numbers, use DOI for citation !!

11