|  Role of Entropy in Determining the Phase Behavior of Protein Solutions
Induced by Multivalent Ions†
---|---
|  Anil Kumar Sahoo,∗abc Frank Schreiber,d Roland R. Netz,ce and Prabal K.
Maiti∗a
|  Recent experiments have reported lower critical solution temperature (LCST)
phase behavior of aqueous solutions of proteins induced by multivalent ions,
where the solution phase separates upon heating. This phenomenon is linked to
complex hydration effects that result in a net entropy gain upon phase
separation. To decipher the underlying molecular mechanism, we use all-atom
molecular dynamics simulations along with the two-phase thermodynamic method
for entropy calculation. Based on simulations of a single BSA protein in
various salt solutions (NaCl, CaCl2, MgCl2, and YCl3) at temperatures ($T$)
ranging 283–323 K, we find that the cation–protein binding affinity increases
with $T$, reflecting its thermodynamic driving force to be entropic in origin.
We show that in the cation binding process, many tightly bound water molecules
from the solvation shells of a cation and the protein are released to the
bulk, resulting in entropy gain. To rationalize the LCST behavior, we
calculate the $\zeta$-potential that shows charge inversion of the protein for
solutions containing multivalent ions. The $\zeta$-potential increases with
$T$. Performing simulations of two BSA proteins, we demonstrate that the
protein–protein binding is mediated by multiple cation bridges and involves
similar dehydration effects that cause a large entropy gain which more than
compensates for rotational and translational entropy losses of the proteins.
Thus, the LCST behavior is entropy-driven, but the associated solvation
effects are markedly different from hydrophobic hydration. Our findings have
direct implications for tuning the phase behavior of biological and soft-
matter systems, e.g., protein condensation and crystallization.

##

††footnotetext: a Center for Condensed Matter Theory, Department of Physics,
Indian Institute of Science, Bangalore-560012, India. E-mail:
<EMAIL_ADDRESS><EMAIL_ADDRESS>b Max Planck
Institute of Colloids and Interfaces, Am Mühlenberg 1, 14476 Potsdam,
Germany.††footnotetext: c Fachbereich Physik, Freie Universität Berlin,
Arnimallee 14, 14195 Berlin, Germany.††footnotetext: d Institute for Applied
Physics, University of Tübingen, 72076 Tübingen, Germany.††footnotetext: e
Department of Physics, Indian Institute of Science, Bangalore-560012,
India.††footnotetext: † Electronic Supplementary Information (ESI) available:
[Distribution of anions around the protein; electrostatic potential map and
cation/anion number density map; cation binding/unbinding kinetics; time
series of the number of protein-bound water molecules; radial distribution
functions for oxygen of water around cations; entropy contribution of water in
the second solvation shell of cations; comparison of entropy obtained from the
2PT method and the temperature derivative of free energy; the $T$-dependence
of binding energy; $T$ and salt concentration dependence of the dielectic
constant; the $T$-dependence of entropy and energy of ion binding excluding
the contribution due to dehydration of the second solvation shell of the ion;
note on the $T$-dependence of free energy; theory of the 2PT method; note on
volume calculation; derivation of the expression for $\zeta$-potential (Eq.
9); table for comparison between computed and experimental values for
structural parameters and entropy of ion hydration; and table for water
released from protein and cation surfaces in a binding process]. See DOI:
10.1039/cXsm00000x/

## Introduction

Ions play an important role in many biophysical processes, e.g., allosteric
regulation, enzymatic activity, DNA condensation, and protein solubility and
crystallization. Starting from the pioneering works by Hofmeister, there has
been immense progress made to better understand ion–protein interactions.1, 2
In recent years, due to various applications in biology, medicine and physics,
there is increasing interest to tune and control the phase behavior of protein
solutions using multivalent ions.3 Diverse phenomena induced by multivalent
ions have been realized in experiments. These include: (i) reentrant
condensation of proteins in bulk solution4 as well as reentrant surface-
adsorption of proteins5 by varying the concentration of Y3+ or other trivalent
cations, (ii) pathway-controlled protein crystallization,6 (iii) clustering,7
(iv) liquid–liquid phase separation,7, 8 and (v) lower critical solution
temperature (LCST) phase behavior.9 Although many aspects regarding
ion–protein interactions have been qualitatively understood, a fundamental and
quantitative understanding is required for further developments in this field.

Of particular interest is the LCST phase behavior for a solution of bovine
serum albumin (BSA) proteins in the presence of Y3+ ions.9 At low
temperatures, the proteins remain well dispersed in solution, whereas upon
increasing temperature up to 300 K, the proteins attract each other, and the
solution separates into protein-rich and protein-poor phases. We note that
aggregation of proteins can also be caused by thermal denaturation, but in the
experiments Matsarskaia et al.9 stayed well below the protein denaturation
temperature and observed LCST behavior only for solutions containing trivalent
ions.10 This precludes denaturation as a mechanism and suggests that the LCST
behavior is related to ion-mediated protein aggregation.

It has been suggested that the LCST behavior is due to the combination of
effects associated with the solvation of the protein and the multivalent ions,
and that entropy is the driving force.9 However, the molecular mechanism of
the LCST behavior has not been quantitatively identified. A quantitative
understanding of the thermodynamics of this process requires an accurate
estimation of various entropy contributions associated with the ion–protein
complex formation and the subsequent ion-mediated protein–protein aggregation.
The total entropy change includes entropy costs due to (i) hindrance in the
translation of a protein-bound ion, (ii) restrictions on the translational and
rotational motions of proteins, (iii) hydration/dehydration of the protein and
ions, and (iv) conformational changes of the protein. The latter is mainly
important for metalloregulatory allosteric proteins. Quantifying all these
entropy contributions in experiments remains a daunting task, even with the
present-day techniques that provide residue-level dynamic information.11 In
this regard, molecular simulations12, 13 along with accurate and robust
entropy calculation techniques provide an alternative and reliable approach.

To understand the mechanistic details and the thermodynamic driving force for
the intriguing phenomena related to ion-mediated protein–protein interactions,
we have performed large-scale molecular dynamics (MD) simulations of a single
and two BSA proteins in chloride salt solutions of Y3+ and several other
cations found in physiological conditions, such as Na+, Ca2+, and Mg2+ in the
temperature range of 283–323 K. The simulation details are presented in the
Methods section. A snapshot of the initial configuration of the simulated
single-protein system is shown in Figure 1A. We investigate the specific
nature of ion–protein interactions and quantify the free energy, various
entropy contributions as well as electrostatics of the system. Our study
reveals crucial solvation/desolvation phenomena giving rise to an entropic
driving force for ion–protein binding, in contrast to common expectations.
From simulations of the systems involving two BSA proteins, it is found that
Y3+ ions link the two proteins to form a dimer. Hence, the process of ion-
mediated protein-protein binding is argued to be entropy-driven, as a large
number of tightly bound water molecules are released from the proteins and the
mediating cations’ surfaces to the bulk solution.

## Results

### BSA Protein–Ion Interaction and Ion Binding Kinetics

To investigate the nature of ion–protein interactions, we calculate the number
distribution of ions $N(r)$ along the protein’s surface-normal direction.
$N(r)$ for the cations are shown in Figure 1B, while $N(r)$ for Cl- ion for
the different ionic solutions are plotted in Figure S1 in the supporting
information (SI). We find that cations are mostly present near the protein,
with the relative propensity of binding showing the following trend:
monovalent < divalent < trivalent. These cations predominantly pair with the
negatively charged carboxylate groups of aspartate and glutamate surface
residues of the protein. Interestingly, even in NaCl solution, Cl- ions are
found to be largely present near the protein, and the number of Cl- ions
present near the protein decreases in the following order: YCl3 > MgCl2
$\approx$ CaCl2 > NaCl (Figure S1 in the SI). This suggests that Cl- ions
interact with the –NH${}_{3}^{+}$ group of the protein surface residues, and
also interact, via ion-pair formation, with the cations present in the
vicinity of the protein.

A protein surface is, however, far from uniform and if some extended patches
are present on its surface, strong affinity of multivalent ions is expected
even if the net charge of the protein is small or even of the opposite sign.14
We indeed find a positively charged patch and a few extended negatively
charged patches from the electrostatic potential map for BSA (Figure S2A in
the SI). We find higher density of cations (anions) near negatively
(positively) charged patches even for monovalent ions (Figure S2B,C in the
SI).

Fig. 1: (A) Snapshot for the starting simulation box (of size
13.2$\times$13.2$\times$13.1 nm3) containing a single BSA protein in 30 mM
YCl3 solution. The protein (orange) is shown in cartoon representation. Y3+
(blue) and Cl- (red) ions are represented as VDW spheres. Water molecules are
represented as a continuum in semi-transparent mode for clarity. Note that the
system was also simulated in other salt solutions (MgCl2, CaCl2, and NaCl) of
the same ionic strength as in the case of YCl3. (B) For each type of cation,
the total number of cations $N$ found within shells of width $dr$ $=0.1$ Å,
present at a shortest distance $r$ from the protein surface, is plotted as a
function of $r$ for the simulation performed at 303 K. $N(r)$ for each cation
type is averaged over the last 100 ns of the simulation time. Note that $N(r)$
is enhanced near the protein surface, representing the strong affinity of
cations for the negatively charged protein ($-16$ $e$ at pH 7).

To check how tightly the cations are bound to the protein, we monitor their
binding/unbinding kinetics. An ion is defined as bound if it is within a
cutoff distance $r_{c}$ from any atom of the protein, otherwise the ion is
unbound or free. From the $N(r)$ plot in Figure 1B, $r_{c}$${}^{\prime}s$ for
the different cations are chosen as 2.8 Å (Na+), 2.7 Å (Ca2+), 2.3 Å (Mg2+),
and 2.5 Å (Y3+). We find intermittent binding/unbinding events for both Na+
and Ca2+ ions (Figure S3 in the SI). While the binding/unbinding events for
Na+ ions are frequent, prolonged bindings are observed for Ca2+ ions. For
these two cation types, the binding time, i.e., the duration for which an ion
remains bound once it comes within distance of $r_{c}$ from the protein, is
broadly distributed, owing to the surface heterogeneity of the protein. In
contrast, only one unbinding event is observed for Mg2+ within 1.27 $\mu$s,
whereas no unbinding of Y3+ is seen within 1.45 $\mu$s (Figure S3 in the SI).
As the water escape time in the first solvation shell of Mg2+ is $\sim$1.5
$\mu$s,15 it presumably requires very long simulations (100 $\mu$s to 10 ms)
to obtain sufficient unbinding statistics for Mg2+ and Y3+ ions. Performing
such long, all-atom simulations for our system is out of reach of our
computational capabilities.

Fig. 2: Time series of the total number of protein-bound cations ($N_{bI}$) at
several temperatures for Na+ (A), Ca2+ (B), Mg2+ (C) and Y3+ (D) ions. For the
multivalent cations, $N_{bI}$ increases upon increasing temperature.

For each cation type, the total number of protein-bound cations, $N_{bI}$, is
plotted as a function of the simulation time at three different temperatures
in Figure 2. No ion is bound to the protein at the beginning of a simulation,
and $N_{bI}$ gradually increases with the simulation time. $N_{bI}$ eventually
reaches a saturation value, at a time required for equilibration of the ion
distribution around the protein. This ion equilibration time differs for each
cation, which can be rationalized by considering the ion–water exchange
kinetics that strongly depends on the cation’s charge and size.15
Counterintuitively, we find from Figure 2 that $N_{bI}$ increases with
increasing temperature. This effect is prominent for all the cations, except
Na+. In contrast, the number of protein-bound water, i.e., the total number of
water molecules present within 3 Å from the protein surface decreases with the
increase in temperature as expected (Figure S4 in the SI). Although an
increase in the binding affinity of any two objects by raising the temperature
is not new—e.g., hydrophobic interaction strength increases with
temperature,16 it is surprising to be observed in a system involving strong
electrostatic interactions and can be rationalized by the temperature
dependence of dielectric and hydration effects.17, 18 For a quantitative
understanding of this, we calculate various thermodynamic quantities such as
the free energy, enthalpy, and various entropy contributions as discussed
below.

### Thermodynamics of Cation Binding to the Protein

The free energy of a cation binding to the protein, $\Delta G_{b}$, for
temperatures in the range of 283–323 K is shown in Figure 3A (see Methods for
the calculation details). For each cation type, $\Delta G_{b}$ is always
negative, and its magnitude increases with the increase in temperature.
$|\Delta G_{b}|$ follows the trend: Na+ < Ca2+ $\approx$ Mg2+ < Y3+. By
changing temperature from 283 K to 323 K, we see the highest change in $\Delta
G_{b}$ for Y3+ ($-1.21$ kcal/mol), whereas the least change is observed for
Na+ binding ($-0.52$ kcal/mol). The changes in $\Delta G_{b}$ for Ca2+ and
Mg2+ ions are $-0.71$ and $-1.03$ kcal/mol, respectively.

The increase in binding affinity of the cations with solely increasing
temperature (Figure 3A) cannot be explained by considering the energy of
binding, for purely thermodynamic reasons, as described in the SI, section 1.
Further, it should be noted that since the dielectric constant of water
decreases as $\sim T^{-3/2}$, any electrostatic interaction in water is
predominantly entropic in nature.17, 18 Therefore, entropy must be playing a
dominant role here.

Fig. 3: Temperature dependence of (A) the free energy, $\Delta G_{b}$, and (B)
the total entropy contribution, $-T\Delta S_{b}$, for each cation binding to
the protein. Error bars represent the standard deviation. The different lines
are for guiding the eye. The experimental binding free energies for Y3+ at
different temperatures shown in (A) are taken from Ref.9

The binding free energy for an ion is given by

$\Delta G_{b}(T)=\Delta E_{b}(T)-T\Delta S_{b}(T),$ (1)

where $\Delta E_{b}$ and $\Delta S_{b}$ are the energy and entropy of binding,
respectively and $T$ is the temperature. For the calculation of $\Delta
S_{b}$, one needs to correctly account for “hydration effects” associated with
the ion binding process, such as partial desolvation of both the protein and
ion. The radial distribution functions for water molecules around a cation,
both free in solution and bound to the protein surface, clearly show partial
dehydration of the first and second solvation shells (SS${}^{\prime}s$) of
each cation (Figure S5 in the SI). $\Delta S_{b}$ in Eq. 1 consists of three
terms—the loss in entropy of a protein-bound ion ($\Delta S_{P,I}$), the gain
in entropy due to release of tightly-bound water molecules from the first and
second SS${}^{\prime}s$ of the ion ($\Delta S_{I,W}$), and the gain in entropy
of water molecules released to the bulk due to desolvation of the protein
surface residue where the ion binds ($\Delta S_{P,W}$). Together, it can be
written as

$\Delta S_{b}=\Delta S_{P,I}-\Delta S_{I,W}-\Delta S_{P,W}.$ (2)

We have used the two-phase thermodynamic (2PT) method19, 20, 21 to calculate
all the terms on the right hand side of Eq. 2. The theory of the 2PT method is
described in the SI, section 2 and calculation details are given in the
Methods section. We first validate the 2PT method for ionic solutions by
reproducing from the simulation data the experimental ion hydration entropy
($\Delta S_{hyd}$) for the different ion types (see Table S1 in the SI). Then,
we proceed, using 2PT, with calculations of the entropy differences for a
protein-bound ion, a protein-bound water, and a water in the first SS of the
cation as shown in Figure 4, as well as the entropy difference of a water in
the second SS of the cation as shown in Figure S6 in the SI. Note that for
calculations of the various entropy contributions shown in Figures 4 and S6,
the reference values are taken as the respective absolute entropies in the
bulk water. $\Delta S_{I,W}$ in Eq. 2 is then calculated by multiplying the
per water entropy differences with the corresponding numbers of water
molecules released in the partial dehydration of both the first and second
SS${}^{\prime}s$ of the cation (values given in Table S2 in the SI), and
adding both terms. Similarly, $\Delta S_{P,W}$ in Eq. 2 is evaluated by
multiplying the per water entropy difference with the number of water
molecules released in the partial dehydration of the protein surface residue
(values given in Table S2 in the SI). From Figure 4, we see that the entropy
loss of a protein-bound cation is more than compensated by the entropy gain of
water molecules released to bulk by the partial dehydration of both the cation
and protein. The cation desolvation entropy contributes the highest to the
thermodynamics of protein–ion binding for all the multivalent cations, whereas
both the protein and ion desolvation entropies contribute equally for Na+
binding.

The total entropy contribution ($-T\Delta S_{b}$ in Eq. 1) for each cation
plotted in Figure 3B is always negative, and it decreases (becomes more
negative) with increasing temperature. We have also estimated $-T\Delta S_{b}$
from the temperature dependence of $\Delta G_{b}$ using the thermodynamic
relation $\Delta S_{b}=-\partial\Delta G_{b}/\partial T$, and we find that the
temperature dependence trend is the same as obtained from the 2PT method,
though the values obtained from both methods match only semi-quantitatively
(see Figure S7 in the SI for comparison). For each cation, $-T\Delta S_{b}$ is
more negative than the binding free energy $\Delta G_{b}$ throughout the
temperature range studied in this work (Figure 3). Therefore, the process of a
cation binding to the protein is entropy driven. The above observations, in
particular, explain the enhancement of the protein-binding affinity of a
multivalent cation with increasing temperature (Figure 2). The total entropy
contribution as shown in Figure 3B is the highest for Y3+ across the whole
temperature range, followed by Ca2+ > Mg2+ $\approx$ Na+—representing a
delicate dependency of entropy on the charge and size of a cation. Note that
though the entropy contribution due to a water molecule released from Mg2+ is
more than that for Na+ and Ca2+ ions (Figure 4), the altered trend in $\Delta
S_{b}$ for Mg2+ in Figure 3B is rationalized by the lower number of water
molecules released in the process of a Mg2+ ion binding, compared to that for
Na+ and Ca2+ bindings.

Fig. 4: Entropy contributions $\Delta S_{P,I}$ (red), $\Delta S_{I,W}$
(black), $\Delta S_{P,W}$ (blue) from Eq. 2 as a function of temperature for
the protein in NaCl (A), CaCl2 (B), MgCl2 (C), and YCl3 (D) solutions. Error
bars represent the standard deviation. The lines are for guiding the eye.

The large (and negative) value of the entropy contribution, $-T\Delta S_{b}$,
must be partially compensated by a positive binding energy $\Delta E_{b}$ to
result in a small (and negative) value of the binding free energy $\Delta
G_{b}$. $\Delta E_{b}$, calculated by using the thermodynamic relation given
in Eq. 1, is plotted as a function of temperature in Figure S8 in the SI.
$\Delta E_{b}$ is positive throughout the temperature range, in agreement with
the experiment,9 but is comparable to the magnitude of $-T\Delta S_{b}$. The
increase in $\Delta E_{b}$ with temperature (Figure S8) can be rationalized by
the enhancement in the electrostatic interaction strength due to the decrease
in the water dielectric constant $\epsilon$, as explained below.17 The
electrostatic free energy $\Delta G\propto\epsilon^{-1}$ and $\epsilon\propto
T^{-\alpha}$, thus $\Delta G=-CT^{\alpha}$. Here $C$ is a constant and the
negative sign is due to $\Delta G<0$ in our case. The entropy follows as
$\Delta S=-\partial\Delta G/\partial T=C\alpha T^{\alpha-1}$. The internal
energy results as $\Delta E=\Delta G+T\Delta S=-CT^{\alpha}+C\alpha
T^{\alpha}=C(\alpha-1)T^{\alpha}$. As long as the exponent $\alpha>1$, $\Delta
E$ is always positive and increases as $T^{\alpha}$. For pure water $\alpha>1$
at all temperatures (Figure S9 in the SI). Although $\alpha$ slightly
decreases with the addition of salt (viz. 1 M NaCl solution in Figure S9),
$\alpha$ is significantly greater than 1 for the temperature regime
investigated in our simulations, which explains the observed increase in
$\Delta E_{b}$ with increasing temperature.

The temperature-dependent increase in $\Delta E_{b}$ follows the trend: Y3+ >
Ca2+ > Na+ $\approx$ Mg2+ (Figure S8 in the SI). By changing temperature from
283 K to 323 K, the change in $\Delta E_{b}$ for Na+, Ca2+, Mg2+, and Y3+ is
found to be 1.82, 2.86, 1.71, and 5.11 kcal/mol, respectively.

The large value of $\Delta E_{b}$ can be understood by considering the
energetic penalties associated with the desolvation of both the protein and
cation. For example, $\Delta E_{b}$ for Y3+ ion at 300 K decreases from 7.50
to 3.71 kcal/mol if we exclude the contribution due to the dehydration of the
second SS of Y3+ (Figure S10B in the SI). Figure S10 also highlights that the
effect of the second SS is significant for the accurate description of
solvation thermodynamics of cations, and cannot be neglected even for
monovalent ions, e.g., Na+.

### Preferential Interaction Coefficients

The interaction of ions with proteins, whether these are enriched or depleted
from the protein surface, can be quantified by experimentally measuring the
preferential interaction coefficient $\Gamma_{23}$. The thermodynamic
definition of $\Gamma_{23}$ is the change in chemical potential of the protein
due to the addition of ions; it can also be expressed as the change in ion
concentration to maintain constant chemical potential when a protein is added
to the solution:22

$\Gamma_{23}=-\left(\frac{\partial\mu_{2}}{\partial\mu_{3}}\right)_{m_{2},T,P}=-\left(\frac{\partial
m_{3}}{\partial m_{2}}\right)_{\mu_{3},T,P},$ (3)

where $\mu$ is the chemical potential, $m$ is the molal concentration, and the
subscripts 1, 2, and 3 stand for water, protein, and ion, respectively. Record
et al.,23 based on the molal concentration definition, developed a two-domain
molecular model for the estimation of $\Gamma_{23}$ in terms of the difference
in ion concentration in the local domain near the protein surface and the bulk
solution as follows:

$\Gamma_{23}=\bigg{\langle}N_{3}^{local}-N_{1}^{local}\bigg{[}\frac{N_{3}^{bulk}}{N_{1}^{bulk}}\bigg{]}\bigg{\rangle},$
(4)

where $N_{i}$ is the number of molecules of type $i$ and $\langle\cdot\rangle$
represents the time average. For the calculation of $\Gamma_{23}$ using Eq. 4
a boundary or a distance cutoff needs to be chosen for defining the local and
bulk domain, but the choice is arbitrary. $\Gamma_{23}$ is instead estimated
at each value of $r$, the distance from the protein surface, asuming that it
is the boundary: $\Gamma_{23}(r)=\langle
N_{3}(r)-N_{1}(r)[N_{3}^{bulk}/N_{1}^{bulk}]\rangle$. The distance $r^{*}$
after which $\Gamma_{23}(r)$ becomes constant is defined as the actual
boundary. In our simulations the total numbers of water molecules ($N_{1}$)
and ions ($N_{3}$) are constant, thus the above expression for
$\Gamma_{23}(r)$ is further simplified as:24

$\Gamma_{23}(r)=\bigg{\langle}N_{3}(r)-N_{1}(r)\bigg{[}\frac{N_{3}-N_{3}(r)}{N_{1}-N_{1}(r)}\bigg{]}\bigg{\rangle}.$
(5)

In a salt solution, cations and anions are distributed around the protein. We
obtain preferential interaction parameters for the cation $\Gamma_{2,+3}(r)$
and anion $\Gamma_{2,-3}(r)$ separately by using $N_{+/-3}(r)$ as the cation
or anion distribution, respectively in Eq. 5. $\Gamma_{2,+3}(r)$ and
$\Gamma_{2,-3}(r)$ are shown for different salt solutions in Figure 5.
Experimentally, it is impossible, however, to separate the cationic and
anionic contribution to the measured value of $\Gamma_{23}$ for a salt
solution. For a salt of monovalent cation and anion, the preferential
interaction parameter is given by23

$\Gamma_{23}=\frac{1}{2}\big{(}\Gamma_{2,+3}+\Gamma_{2,-3}-|Q_{2}|\big{)},$
(6)

where $|Q_{2}|$ is the protein’s net charge that is subtracted from
$\Gamma_{23}$, as $Q_{2}$ counterions (cations in case of BSA protein) are
accumulated near the protein surface to neutralize its charge and do not
contribute to the preferential interaction. For a salt of multivalent
cation/anion, it is straight-forward to generalize Eq. 6 by scaling
$\Gamma_{2,+3}(r)$, $\Gamma_{2,-3}(r)$, and $Q_{2}$ with valency of the anion
$z_{-}$, valency of the cation $z_{+}$, and charge on the counterion $z_{+}$,
respectively.

$\Gamma_{23}=\frac{1}{2}\bigg{(}\frac{\Gamma_{2,+3}}{z_{-}}+\frac{\Gamma_{2,-3}}{z_{+}}-\frac{|Q_{2}|}{z_{+}}\bigg{)}.$
(7)

For the BSA protein, using in Eq. 7 $\Gamma_{2,+3}$ and $\Gamma_{2,-3}$ at
$r^{*}=17$ Å (by which all curves reach their respective saturation values as
seen in Figure 5), we obtain preferential interaction coefficients for
different salts: NaCl ($\Gamma_{23}=2.44$), CaCl2 ($\Gamma_{23}=15.67$), MgCl2
($\Gamma_{23}=19.83$), and YCl3 ($\Gamma_{23}=26.87$). Positive values of
$\Gamma_{23}$ for all the different salt types reflect that these salt ions
are attracted towards the protein surface. For salt containing multivalent
ions, $\Gamma_{23}$ is significantly larger than that for NaCl, which suggests
that addition of trivalent ions in the protein solution affects the solution
stability25 and stabilizes protein dimer formation as seen in our simulations.

Fig. 5: Preferential interaction coefficients for cations $\Gamma_{2,+3}$ and
anions $\Gamma_{2,-3}$ of the different salt solutions at 303 K as a function
of the boundary distance cutoff $r$. $\Gamma_{2,+3}(r)$ reaches a saturation
value rather quickly compared to $\Gamma_{2,-3}(r)$, for solutions containing
multivalent cations.

###  $\zeta$-Potential of the Protein and the Protein–Protein Interaction

$\zeta$-potential measurements for a protein in an ionic solution report on
charge compensation by the counterions and thus have direct implications for
protein–protein association and the phase behavior of the solution.
$\zeta$-potentials are defined by the electrophoretic mobility.10, 9, 26 From
the simulation data, we calculate the surface potential at one ionic diameter
away from the protein surface (see Methods). Note that the surface potential
typically serves as a good approximation for the $\zeta$-potential for
proteins and colloidal systems; however, the surface and $\zeta$-potential
values might differ significantly for extended surfaces with high surface
charge densities.27

As shown in Figure 6, the $\zeta$-potential of the protein in the NaCl
solution is negative at all temperatures, as expected based on the protein net
charge of $-16$ $e$. In contrast, the $\zeta$-potential is positive for all
multivalent cation-chloride solutions at all temperatures, indicating sign
reversal of the effective charge of the protein (Figure 6). This charge
inversion phenomenon in the presence of multivalent cations is due to strong
interactions of the cations predominantly with the COO- groups of the
protein’s surface residues and can be rationalized by considering strong
charge-charge correlations.28 Note that similar charge reversal of proteins in
the presence of trivalent cations has also been reported both in experiments4,
26, 9 and simulations13 as well as in a coarse-grained analytical model.29 As
shown in Figure 6, with the increase in temperature, the $\zeta$-potential of
the protein increases for all cation types; the highest change is seen for
MgCl2, whereas the effect is minimal for NaCl. The $\zeta$-potential of the
protein at 283 K is higher in CaCl2 solution than in MgCl2 solution, and vice
versa at 323 K. These observations are consistent with the trends for the
temperature dependence of binding free energies of the different cations
(Figure 3A).

Fig. 6: Temperature dependence of the $\zeta$-potential of the protein in
various salt solutions. The lines are guide to the eye only.

### Protein–Protein Binding Mediated by Cation Bridges

Protein aggregation seen in experiments9 was hypothesized to be mediated by
cation bridges.6 To explicitly demonstrate the multivalent ion-mediated
protein–protein binding, we have performed three independent simulations with
different orientations of two BSA proteins in YCl3 solution, as shown in
Figure 7A (left panel). In every simulation, we find that two BSA proteins,
which are initially placed far apart, approach each other (see the timeseries
of the total number of inter-protein residue–residue contacts in Figure 7B)
and eventually form a dimer mediated by 1–5 Y3+ ions (see snapshots in Figure
7A [middle panel]). The Y3+ ion bridges remain stable over a 1 $\mu$s
timescale, as evident from the time series plot for the number of bridging
cations (Figure 7C). A Y3+ ion bridge is stabilized by coordination of
multiple caboxylate groups of each protein with the cation, as evident from
snapshots in Figure 7A (right panel). Note that even for the stable, Y3+ ion-
bridged protein dimer complex, the relative orientation between the two
proteins changes over time but very slowly (see the orientational
autocorrelation function in Figure 7D). This reveals the conformational
flexibility of the protein dimer complex.

To compare monovalent and multivalent cations, we have also simulated the
above three systems (shown in Figure 7A) in NaCl solution, at the equivalent
ionic strength as for YCl3 solution. In sharp contrast to the case of YCl3, we
find that Na+ ion bridges between two BSA proteins form transiently and remain
stable only for 1–20 ns (Figure 7E). These results demonstrate the need for
multivalent ions in protein cluster formations, in agreement with
experiments.4

Fig. 7: Results for two BSA proteins (A–D) in 30 mM YCl3 and (E) 180 mM Nacl
solution at 303 K. (A) The initial and final configuration of two BSA proteins
(represented in orange and tan). Bridging Y3+ ions in the final structure are
shown as blue spheres, and amino acid residues of the two proteins involved in
the formation of ion bridges are highlighted in the ball–stick representation
(see the right panel for the zoomed-in version). Water and other ions are
omitted for clarity. (B) Time series of the total number of inter protein
residue–residue contacts. (C) Time series of the total number of Y3+ ion
bridges that link the two proteins to form a dimer. (D) The average
autocorrelation function for the relative orientation angles between the
principal axes of two proteins as a function of time delay. (E) Time series of
the total number of Na+ ion bridges that form transiently between two BSA
proteins in 180 mM NaCl solution at 303 K. All observables are defined in
Methods.

## Discussion

The temperature behavior of the $\zeta$-potential found in Figure 6 is in
qualitative agreement with the experiments.9 As the $\zeta$-potential is
influenced by the number of surface-bound ions and the binding affinity of
ions increases with temperature, the $\zeta$-potential is expected to increase
with temperature irrespective of the salt concentration of the solution.
$\zeta$-potential values estimated from our simulations, however, are larger
than that reported in the experiments9 presumably because of the YCl3
concentration difference. The $\zeta$-potential increases with increasing YCl3
concentration as found in experiments,9 thus we expect the simulation and
experimental results to match if the same salt concentration is used.

It should be noted that the YCl3 concentration used in our simulations is 30
mM which is much higher than the 1 mM concentration used in the experiment. A
direct comparison between all-atom simulations and the experiments9 at low
concentration of multivalent ions is rather difficult for the following
reason. We consider a higher YCl3 ion concentration in our simulations in
order to obtain statistically converged results with sufficient number of
ions. Obtaining well-converged results for proteins at low salt concentrations
with enough number of ions would require significantly larger system sizes.
Simulating such large systems is very demanding at the all-atom level, but it
is feasible at a coarse-grained level as shown in a recent study.13 However,
the solvation effects, which are crucial for the accurate prediction of
protein–ion binding thermodynamics, are not properly taken into account in
such coarse-grained simulations.

The LCST phase behavior found in experiments9 can be rationalized by the
temperature dependence of the $\zeta$-potential and the microscopic picture
emerges from our simulations. For sufficiently low Y3+ concentration, at low
temperatures due to the reduced binding affinity of counterions, the
$\zeta$-potential is expected to become negative (and large) and the proteins
are expected to repel each other, keeping the solution stable. With increasing
temperature, counterion binding affinity for the protein increases, and hence
the $\zeta$-potential increases and becomes positive at a sufficiently high
temperature. In the temperature range (293–313 K) where the $\zeta$-potential
is small ($-$5 to $+$5 mV),9 the proteins are predicted to attract each other,
eventually causing the solution to phase separate into protein-rich and
protein-poor phases.

The protein–protein binding at a low concentration of multivalent ions occurs
via cation bridging, as shown in Figure 7, as well as suggested from
experiments.6, 9 A cation bridge formation—like the first step of a cation
binding to the protein—requires desolvation of both the protein-bound cation
and the surface residue of another protein that will bind to the cation. These
processes involve the release of many tightly-bound water molecules to the
bulk that results in a significant entropy gain, which contributes at least
10–15 kcal/mol (depending on the temperature) to the total free energy, as
shown in Figure 3B for a Y3+ ion binding. As multiple cation bridges are
formed in a protein–protein binding (see Figure 7, and also found in
experiments6), the net entropy gain due to cation and protein desolvation more
than compensates for the translational and rotational entropy losses of the
proteins during protein–protein binding. Therefore, the LCST phase behavior9
is entropy-driven.

## Conclusions

In summary, by performing fully atomistic MD simulations of a BSA protein in
different cation-chloride solutions (NaCl, CaCl2, MgCl2, and YCl3) and by
calculating various entropy contributions, we demonstrate that multivalent
cation binding to the protein is an entropy-driven phenomenon. The loss in
entropy of a protein-bound cation is more than compensated by the entropy gain
of water molecules due to the partial dehydration of both the cation and the
cation-bound surface residue of the protein. A particularly interesting
observation is the significant difference in the binding/unbinding kinetics of
Ca2+ and Mg2+ ions (see Figure S3)—although having comparable binding free
energies (Figure 3A), which can be related to the recent finding that the
ion–water exchange kinetics strongly depends on the size of a cation.15 It
will thus be interesting to investigate in future simulation studies the
universality of the ion size dependence of ion–protein binding kinetics and
thermodynamics.

The $\zeta$-potential calculation shows charge inversion of the protein in all
solutions containing multivalent cations, but not in the monovalent NaCl
solution (Figure 6). The LCST phase behavior observed in the experiment9 can
be rationalized by considering the temperature-dependent increase in the
$\zeta$-potential of the protein and the associated charge inversion
phenomenon. The protein–protein interaction involves: (i) the ion binding to
the protein, and (ii) the subsequent protein–protein binding by cation
bridging (Figure 7). In both processes many tightly-bound water molecules are
released to the bulk, which results in a thermodynamic driving force for the
LCST behavior that is entropic in nature, in agreement with the experiment.9

This work shows that similarly to hydrophobic association, entropy plays a
pivotal role in systems involving strong electrostatic interactions, revealing
intriguing hydration and dielectric effects. Our results are important for the
basic understanding of ion effects in soft matter and biology, and the
insights gained here will be useful in studies of ion-mediated surface
adsorption and crystallization of proteins. Moreover, molecular-level
understanding of interactions of heavy metals—usually not found in healthy
cells—with different biomolecules, as studied here, can provide insights for
carcinogenicity and neurotoxicity induced by exposure to such environmental
contaminants.

## Methods

### Model Building and Force Field Parameters

The initial structure of BSA protein was obtained from the crystal structure
available in the protein data bank (PDB ID: 3V03). The charge or protonation
state of each residue of the protein was chosen at neutral pH 7 depending on
the residue’s pKa value, and the assigned charges were fixed over the
simulation time. Note, however, that pKa depends on the ionic strength
(through the activity coefficients), and the reported pKa values of amino
acids are typically determined in a solution of high ionic strength.30 In
particular, the apparent pKa values of carboxyl groups shift up slightly in
the presence of multivalent cations at low salt concentrations. If the pH of
the solution differs from 7 in an experiment due to the CO2 content in air
(which lowers the water pH down to 5.6 31), this could make some of the acidic
peptide groups less charged. But, the experiments described in Matsarskaia et
al.9 were performed in air and in ultrapure (MilliQ, 18.2 Mega Ohm) water
which had previously been degassed under vacuum to eliminate the CO2
contributions. Also, it is known from experiments that the addition of
multivalent metal cations such as Al3+ and Fe3+ changes the pH of the solution
due to hydrolysis of these cations, which can change the charge states of the
protein surface residues. However, this effect is less significant for Y3+
ion.10 Therefore, our assigned fixed charges of the protein residues at pH 7
for the different cations (Na+, Ca2+, Mg2+, and Y3+) is assumed to mimic the
experimental conditions sufficiently well.

The ff14SB force field parameters32 were used for the protein. The system was
solvated with TIP3P33 water model using the xleap module of the AMBER17
tools34 in a way such that there exists at least 17 Å solvation shell in
between the solute and simulation box wall. The final unit shell for
simulation is a rectangular box of size 13.2$\times$13.2$\times$13.1 nm3 that
contains $\sim$200,000 atoms (Figure 1A). The system was simulated in four
different salt solutions, namely NaCl, MgCl2, CaCl2 and YCl3. Depending on the
ion type, an appropriate number of counterions were added to ensure the charge
neutrality of the simulation unit shell. To simulate the system at a specified
salt concentration, enough numbers of counterions/coions, estimated from the
mole fraction of counterions/coions and water, were further added to the
system. For YCl3, the system was simulated at 30 mM salt solution. For the
other salts, the system was simulated at the equivalent ionic strength as in
the case of YCl3, e.g., 180 mM for NaCl. Especially for multivalent ions, the
electronic polarization effect contributes significantly to the total
interaction energy of such an ion with another charged object. The recently
developed Li/Merz ion parameters35, 36, 37 with 12-6-4 Lennard-Jones (LJ)-type
nonbonded interaction terms take care of the electronic polarization effect
and have been shown to well reproduce the experimental measurables, such as
the ion–oxygen (of water) distance, the ion–water coordination number, and the
hydration free energy of mono- and multi-valent ions. We have provided in
Table S1 in the SI the structural parameters and entropy of ion hydration for
the different ions calculated from our simulation data, which quantitatively
match with the corresponding experimental values. Therefore, we used Li/Merz
ion parameters for an accurate modeling of the ion–water and ion–protein
interactions.

### MD Simulation Details

All the simulations were performed using the PMEMD module of the AMBER14
package.38 Periodic boundary condition was used for all the simulations. Bonds
involving hydrogen atoms were constrained using the SHAKE algorithm39 that
allowed the use of a time step of 2 fs for the integration of Newton’s
equation of motion. The temperature of the system was maintained using a
Langevin thermostat40 with the collision frequency of 5.0 ps-1. Berendsen weak
coupling method41 was used to apply a pressure of 1 atm with isotropic
position scaling with a pressure relaxation time constant of 2.0 ps. Particle
mesh Ewald42 sum was used to compute long-range electrostatic interactions
with a real space cutoff of 10 Å. van der Waals and direct electrostatic
interactions were truncated at the cutoff. The direct sum non-bonded list was
extended to cutoff + “nonbond skin” (10 + 2 Å).

The solvated systems with harmonic restraints (force constant of 500
kcal/mol/Å2) on the position of each atom of the protein were first subjected
to 2000 steps of steepest descent energy minimization, followed by 1000 steps
of conjugate gradient minimization to remove bad contacts present in the
initially built systems. The restraints on the protein atoms were sequentially
decreased to zero during further 4000 steps of energy minimization. The energy
minimized systems were then slowly heated from 10 K to the desired temperature
in many steps during the first 210 ps of MD simulation. During this time, the
solute particles were restrained to their initial positions using harmonic
restraints with a force constant of 20 kcal/mol/Å2. The first 2 ns of
equilibration simulations were performed in the NPT ensemble to attain the
proper water density. Simulations were then switched to the NVT ensemble for
further production runs of 200–1450 ns, depending on cation types.

### Data Analysis

All the analyses were carried out by using home-written codes and/or the
AMBER17 tools.34 Images were rendered using the Visual Molecular Dynamics
software.43

The free energy of ion binding, $\Delta G_{b}$, was calculated using the
expression

$\Delta G_{b}=-k_{B}T\ln(C_{bI}/C_{fI}),$ (8)

where $k_{B}$ is the Boltzmann constant, and $C_{bI}$ and $C_{fI}$ are the
concentration of bound and free ions, respectively. The expressions for
calculations of the concentrations are $C_{bI}=N_{bI}/V_{s}$ and
$C_{fI}=N_{fI}/V_{f}$, where $V_{s}$ is the volume of the shell around the
protein surface where ions are considered as bound, $N_{fI}$ (= total number
of ions $-N_{bI}$) is the number of free ions, and $V_{f}$ is the free volume
available for ions. The volumes were calculated following the protocol
described in Ref. 44, by using the Gromacs program gmx sasa.45 Further details
on the volume calculation are provided in the SI, section 3. The last 200 ns
data for each ion type (the last 150 ns for Na+) was taken for the calculation
of $\Delta G_{b}$, whereas the rest of the data served for the equilibration.

The reported entropy contributions in Figures 4 and S6 were obtained by
calculating absolute molar entropies for free and protein-bound ions, free and
protein-bound water molecules, and water molecules in the first and second
SS${}^{\prime}s$ of the cation by using the 2PT method.19, 20, 21 To generate
trajectories for 2PT calculations, simulations were restarted after 100–500
ns, depending on the ion type. 3 short (40 ps) NVT trajectories for each
system at each temperature were generated with coordinates and velocities
saved every 4 fs. To calculate the 2PT entropy for bound ions, we performed
the analysis for all bound ions and got the average entropy per ion, similarly
for bound water. In Figures 4 and S6, each point and the corresponding error
bar are the average and standard deviation of 3 different simulations,
respectively.

The surface or $\zeta$-potential was obtained by calculating as a function of
$r$ (the distance from the center of mass of the protein) the electrostatic
potential profile, $\phi(r)$, for the system as follows.46 $\phi(r)$ was
calculated by solving the Poisson equation, i.e, by carrying out a double
integration of the charge density profile, $\rho(r)$, obtained from our MD
simulation by using the following expression.

$\phi(R)-\phi(r)=-\frac{1}{\epsilon}\int_{r}^{R}dr_{1}\frac{1}{r_{1}^{2}}\int_{0}^{r_{1}}dr_{2}r_{2}^{2}\rho(r_{2}),$
(9)

where $R$ (= 65 Å) is the radius of the inscribed sphere within the
rectangular MD simulation box and $\epsilon$ is the dielectric permittivity of
water. A derivation of Eq. 9 is given in the SI, section 4. At a temperature
$T$, $\epsilon$ was calculated from the Bjerrum length, $\lambda_{B}$, of
water (= 7 Å) by using the relation: $\epsilon=e^{2}/4\pi\lambda_{B}k_{B}T$,
where $e$ is the elementary charge. It should, however, be noted that for the
temperature behavior of $\epsilon\propto T^{-\alpha}$, the exponent $\alpha$
in experiments is close to ${}^{3}/_{2}$ but we take it to be 1 which is close
to what is seen in simulations for a rather similar water model.18 So the
above approximation for $\epsilon(T)$ is deemed to be good for our purpose.
Finally, the $\zeta$-potential was obtained as
$\zeta=\phi(R)-\phi(R_{h}+2r_{c})$. Here, the hydrodynamic radius of the
protein, $R_{h}$, was taken to be 36 Å,47 and $r_{c}$${}^{\prime}s$ for the
different cations (Figure 1B) were taken to be 2.8 Å (Na+), 2.7 Å (Ca2+), 2.3
Å (Mg2+) and 2.5 Å (Y3+).

$N_{Contacts}$ shown in Figure 7B is defined as the total number of inter-
protein amino acid residue–residue contacts, and such a contact is counted if
at least one pair of atoms from residues of two different proteins are within
3 Å.

$N_{IonBridges}$ shown in Figure 7C is defined as the total number of ions
bridging two different proteins, and an ion bridge is counted if an ion is
present within 3 Å from both the proteins’ surfaces.

The average orientational autocorrelation function shown in Figure 7D is
defined as

$C(t)=\frac{1}{3}\sum_{i=1}^{3}\frac{\langle\cos\theta_{i}(0)\cdot\cos\theta_{i}(t)\rangle}{\langle\cos\theta_{i}(0)\cdot\cos\theta_{i}(0)\rangle},$
(10)

with $\cos\theta_{i}(t)=\hat{e}_{i}^{A}\cdot\hat{e}_{i}^{B}$, where
$\hat{e}_{i}^{A}$ and $\hat{e}_{i}^{B}$ are the unit vectors along the
principal axes of proteins A and B, respectively, and the angular bracket
represents the average over time origins.

## Conflicts of interest

There are no conflicts to declare.

## Acknowledgements

We thank Profs. Daan Frenkel and Tod Pascal for helpful discussions. We
acknowledge SERC, IISc for the allocation of computing time at the SAHASRAT
machine. A.K.S. thanks MHRD, India for the research fellowship and the Max
Planck Society for financial support via the MaxWater initiative. R.R.N.
thanks Infosys Foundation for support during his stay at IISc, Bangalore.

## Notes and references

  * Jungwirth and Cremer 2014 P. Jungwirth and P. S. Cremer, _Nat. Chem._ , 2014, 6, 261–263
  * Okur _et al._ 2017 H. I. Okur, J. Hladilkova, K. B. Rembert, Y. Cho, J. Heyda, J. Dzubiella, P. S. Cremer and P. Jungwirth, _J. Phys. Chem. B_ , 2017, 121, 1997–2014
  * Matsarskaia _et al._ 2020 O. Matsarskaia, F. Roosen-Runge and F. Schreiber, _ChemPhysChem_ , 2020, 21, 1742
  * Zhang _et al._ 2008 F. Zhang, M. Skoda, R. Jacobs, S. Zorn, R. A. Martin, C. Martin, G. Clark, S. Weggler, A. Hildebrandt, O. Kohlbacher and F. Schreiber, _Phys. Rev. Lett._ , 2008, 101, 148101
  * Fries _et al._ 2017 M. R. Fries, D. Stopper, M. K. Braun, A. Hinderhofer, F. Zhang, R. M. Jacobs, M. W. Skoda, H. Hansen-Goos, R. Roth and F. Schreiber, _Phys. Rev. Lett._ , 2017, 119, 228001
  * Zhang _et al._ 2011 F. Zhang, G. Zocher, A. Sauter, T. Stehle and F. Schreiber, _J. Appl. Crystallogr._ , 2011, 44, 755–762
  * Zhang _et al._ 2012 F. Zhang, F. Roosen-Runge, A. Sauter, R. Roth, M. W. Skoda, R. M. Jacobs, M. Sztucki and F. Schreiber, _Faraday Discuss._ , 2012, 159, 313–325
  * Zhang _et al._ 2012 F. Zhang, R. Roth, M. Wolf, F. Roosen-Runge, M. W. Skoda, R. M. Jacobs, M. Stzucki and F. Schreiber, _Soft Matter_ , 2012, 8, 1313–1316
  * Matsarskaia _et al._ 2016 O. Matsarskaia, M. K. Braun, F. Roosen-Runge, M. Wolf, F. Zhang, R. Roth and F. Schreiber, _J. Phys. Chem. B_ , 2016, 120, 7731–7736
  * Roosen-Runge _et al._ 2013 F. Roosen-Runge, B. S. Heck, F. Zhang, O. Kohlbacher and F. Schreiber, _J. Phys. Chem. B_ , 2013, 117, 5777–5787
  * Capdevila _et al._ 2018 D. A. Capdevila, K. A. Edmonds, G. C. Campanello, H. Wu, G. Gonzalez-Gutierrez and D. P. Giedroc, _J. Am. Chem. Soc._ , 2018, 140, 9108–9119
  * Hess and van der Vegt 2009 B. Hess and N. F. van der Vegt, _Proc. Natl. Acad. Sci. U.S.A._ , 2009, 106, 13296–13300
  * Pasquier _et al._ 2017 C. Pasquier, M. Vazdar, J. Forsman, P. Jungwirth and M. Lund, _J. Phys. Chem. B_ , 2017, 121, 3000–3006
  * Yigit _et al._ 2015 C. Yigit, J. Heyda and J. Dzubiella, _J. Chem. Phys._ , 2015, 143, 08B606_1
  * Lee _et al._ 2017 Y. Lee, D. Thirumalai and C. Hyeon, _J. Am. Chem. Soc._ , 2017, 139, 12334–12337
  * Chandler 2005 D. Chandler, _Nature_ , 2005, 437, 640
  * Israelachvili 2015 J. N. Israelachvili, _Intermolecular and Surface Forces_ , Academic press, 2015
  * Sedlmeier and Netz 2013 F. Sedlmeier and R. R. Netz, _J. Chem. Phys._ , 2013, 138, 03B609
  * Lin _et al._ 2003 S.-T. Lin, M. Blanco and W. A. Goddard III, _J. Chem. Phys._ , 2003, 119, 11792–11805
  * Lin _et al._ 2010 S.-T. Lin, P. K. Maiti and W. A. Goddard III, _J. Phys. Chem. B_ , 2010, 114, 8191–8198
  * Pascal _et al._ 2011 T. A. Pascal, S.-T. Lin and W. A. Goddard III, _Phys. Chem. Chem. Phys._ , 2011, 13, 169–181
  * Pierce _et al._ 2008 V. Pierce, M. Kang, M. Aburi, S. Weerasinghe and P. E. Smith, _Cell Biochem. Biophys._ , 2008, 50, 1–22
  * Record Jr and Anderson 1995 M. T. Record Jr and C. F. Anderson, _Biophys. J._ , 1995, 68, 786–794
  * Shukla _et al._ 2009 D. Shukla, C. Shinde and B. L. Trout, _J. Phys. Chem. B_ , 2009, 113, 12546–12554
  * Zhang _et al._ 2010 F. Zhang, S. Weggler, M. J. Ziller, L. Ianeselli, B. S. Heck, A. Hildebrandt, O. Kohlbacher, M. W. Skoda, R. M. Jacobs and F. Schreiber, _Proteins: Struct., Funct., Bioinf.,_ , 2010, 78, 3450–3457
  * Kubíčková _et al._ 2012 A. Kubíčková, T. Křížek, P. Coufal, M. Vazdar, E. Wernersson, J. Heyda and P. Jungwirth, _Phys. Rev. Lett._ , 2012, 108, 186101
  * Uematsu _et al._ 2018 Y. Uematsu, R. R. Netz and D. J. Bonthuis, _Langmuir_ , 2018, 34, 9097–9113
  * Grosberg _et al._ 2002 A. Y. Grosberg, T. Nguyen and B. Shklovskii, _Rev. Mod. Phys._ , 2002, 74, 329–345
  * Roosen-Runge _et al._ 2014 F. Roosen-Runge, F. Zhang, F. Schreiber and R. Roth, _Sci. Rep._ , 2014, 4, 1–5
  * Reijenga _et al._ 2013 J. Reijenga, A. van Hoof, A. van Loon and B. Teunissen, _Anal. Chem. Insights_ , 2013, 8, ACI–S12304
  * Uematsu _et al._ 2017 Y. Uematsu, D. J. Bonthuis and R. R. Netz, _J. Phys. Chem. Lett._ , 2017, 9, 189–193
  * Maier _et al._ 2015 J. A. Maier, C. Martinez, K. Kasavajhala, L. Wickstrom, K. E. Hauser and C. Simmerling, _J. Chem. Theory Comput._ , 2015, 11, 3696–3713
  * Jorgensen _et al._ 1983 W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey and M. L. Klein, _J. Chem. Phys._ , 1983, 79, 926–935
  * Case _et al._ 2017 D. A. Case, D. Cerutti, I. T.E. Cheatham, T. Darden, R. Duke, T. Giese, H. Gohlke, A. Goetz, D. Greene, N. Homeyer, S. Izadi, A. Kovalenko, T. Lee, S. LeGrand, P. Li, C. Lin, J. Liu, T. Luchko, R. Luo, D. Mermelstein, K. Merz, G. Monard, H. Nguyen, I. Omelyan, A. Onufriev, F. Pan, R. Qi, D. Roe, A. Roitberg, C. Sagui, C. S. W. Botello-Smith, J. Swails, R. Walker, J. Wang, R. Wolf, X. Wu, L. Xiao, D. York and P. Kollman, _AMBER 2017_ , University of California, San Francisco, 2017
  * Li and Merz Jr 2013 P. Li and K. M. Merz Jr, _J. Chem. Theory. Comput._ , 2013, 10, 289–297
  * Li _et al._ 2014 P. Li, L. F. Song and K. M. Merz Jr, _J. Phys. Chem. B_ , 2014, 119, 883–895
  * Li _et al._ 2015 P. Li, L. F. Song and K. M. Merz Jr, _J. Chem. Theory Comput._ , 2015, 11, 1645–1657
  * Case _et al._ 2005 D. A. Case, T. E. Cheatham, T. Darden, H. Gohlke, R. Luo, K. M. Merz, A. Onufriev, C. Simmerling, B. Wang and R. J. Woods, _J. Comput. Chem._ , 2005, 26, 1668–1688
  * Ryckaert _et al._ 1977 J.-P. Ryckaert, G. Ciccotti and H. J. C. Berendsen, _J. Comput. Phys._ , 1977, 23, 327–341
  * Van Gunsteren and Berendsen 1988 W. F. Van Gunsteren and H. J. C. Berendsen, _Mol. Simul._ , 1988, 1, 173–185
  * Berendsen _et al._ 1984 H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, A. DiNola and J. R. Haak, _J. Chem. Phys._ , 1984, 81, 3684–3690
  * Darden _et al._ 1993 T. Darden, D. York and L. Pedersen, _J. Chem. Phys._ , 1993, 98, 10089–10092
  * Humphrey _et al._ 1996 W. Humphrey, A. Dalke and K. Schulten, _J. Mol. Graph._ , 1996, 14, 33–38
  * Becconi _et al._ 2017 O. Becconi, E. Ahlstrand, A. Salis and R. Friedman, _Isr. J. Chem._ , 2017, 57, 403–412
  * Abraham _et al._ 2015 M. J. Abraham, T. Murtola, R. Schulz, S. Páll, J. C. Smith, B. Hess and E. Lindahl, _SoftwareX_ , 2015, 1, 19–25
  * Maiti and Messina 2008 P. K. Maiti and R. Messina, _Macromolecules_ , 2008, 41, 5002–5006
  * Li _et al._ 2012 Y. Li, G. Yang and Z. Mei, _Acta. Pharm. Sin. B_ , 2012, 2, 53–59


# GRADIENT GIBBS MEASURES WITH PERIODIC BOUNDARY LAWS OF A GENERALIZED SOS
MODEL ON A CAYLEY TREE:

F. H. Haydarov, R. A. Ilyasova  F. H. Haydarov1,2,3,4, R. A. Ilyasova4 1\. New
Uzbekistan University, 54, Mustaqillik Ave., Tashkent, 100007, Uzbekistan,
2\. AKFA University, 17, Kichik Halqa Yuli Street, Tashkent city, 100095,
Uzbekistan.  3\. Institute of Mathematics, 9, University str., Tashkent,
100174, Uzbekistan,  4\. National University of Uzbekistan, University str., 4
Olmazor district, Tashkent, 100174, Uzbekistan<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract.

We consider Gradient Gibbs measures corresponding to a periodic boundary law
for a generalized SOS model with spin values from a countable set, on Cayley
trees. On the Cayley tree, detailed information on Gradient Gibbs measures for
models of SOS type are given in [3, 11, 8, 12] and we continue the works for
the generalized SOS model. Namely, in this paper the problem of finding
Gradient Gibbs measures which correspond to periodic boundary laws is reduced
to a functional equation and by solving the equation all Gradient Gibbs
measures with 4 periodic boundary laws are found.

Mathematics Subject Classifications (2010). 60K35 (primary); 82B05, 82B20
(secondary)

Key words. Generalized SOS model, spin values, Cayley tree, gradient Gibbs
measure, periodic boundary law.

##  1\. Preliminaries

The gradient Gibbs measure is a probability measure on the space of gradient
fields defined on a manifold. It is often used in statistical mechanics to
describe the equilibrium states of a system. The gradient Gibbs measure is
derived from the Gibbs measure, which is a probability measure on the space of
field configurations. The critical difference is that the gradient Gibbs
measure focuses on the gradients of the fields rather than the fields
themselves (e.g. [7]).

Specifically, the Gradient Gibbs measure is defined on the set of spin
configurations of a system on a Cayley tree. The Gradient Gibbs measure on a
Cayley tree assigns a probability to each possible spin configuration based on
the energy of that configuration. The energy of a spin configuration is
determined by the interactions between neighboring spins. In the case of a
Cayley tree, each spin is coupled to its nearest neighbors along the edges of
the tree (see [5]).

Mathematically, the gradient Gibbs measure assigns a probability to each
possible configuration of a gradient field on the Cayley tree, based on an
energy function. The energy function typically represents the interactions
between the gradients of a scalar field or a vector field. The probability of
a configuration is proportional to the exponential of the negative energy of
that configuration (e.g. [1, 4, 5, 13, 15]).

The study of random field $\xi_{x}$ from a lattice graph (e.g.,
$\mathbb{Z}^{d}$ or a Cayley tree $\Gamma^{k}$ ) to a measure space
$(E,\mathcal{E})$ is a central component of ergodic theory and statistical
physics. In many classical models from physics (e.g., the Ising model, the
Potts model, the SOS model), $E$ is a finite set (i.e., with a finite
underlying measure $\lambda$ ), and $\xi_{x}$ has a physical interpretation as
the spin of a particle at location $x$ in a crystal lattice (detail in [1, 2,
3, 6, 7, 8, 9, 10, 15, 16]).

Let us give basic definitions and some known facts related to (gradient) Gibbs
measures. The Cayley tree $\Gamma^{k}=(V,L)$ of order $k\geq 1$ is an infinite
tree, i.e. connected and undirected graph without cycles, each vertex of which
has exactly $k+1$ edges. Here $V$ is the set of vertices of $\Gamma^{k}$ va
$L$ is the set of its edges.

Consider models where the spin takes values in the set
$\Phi\subseteq\mathbb{R}_{\infty}^{+}$, and is assigned to the vertices of the
tree. Let $\Omega_{A}=\Phi^{A}$ be the set of all configurations on $A$ and
$\Omega:=\Phi^{V}$. A partial order $\preceq$ on $\Omega$ defined pointwise by
stipulating that $\sigma_{1}\preceq\sigma_{2}$ if and only if
$\sigma_{1}(x)\leq\sigma_{2}(x)$ for all $x\in V$. Thus $(\Omega,\preceq)$ is
a poset, and whenever we consider $\Omega$ as a poset then it will always be
with respect to this partial order. The poset $\Omega$ is complete. Also,
$\Omega$ can be considered as a metric space with respect to the metric
$\rho:\Omega\times\Omega\rightarrow\mathbb{R}^{+}$ given by

$\rho\left(\left\\{\sigma(x_{n})\right\\}_{x_{n}\in
V},\left\\{\sigma^{\prime}(x_{n})\right\\}_{x_{n}\in V}\right)=\sum_{n\geq
0}2^{-n}\mathcal{X}_{\sigma(x_{n})\neq\sigma^{\prime}(x_{n})},$

where $V=\\{x_{0},x_{1},x_{2},....\\}$ and $\mathcal{X}_{A}$ is the indicator
function.

We denote by $\mathcal{N}$ the set of all finite subsets of $V$. For each
$A\in V$ let $\pi_{A}:\Omega\rightarrow\Phi^{A}$ be given by
$\pi_{A}\left(\sigma_{x}\right)_{x\in V}=\left(\sigma_{x}\right)_{x\in A}$ and
let
$\mathcal{C}_{A}=\pi_{A}^{-1}\left(\mathcal{P}\left(\Phi^{A}\right)\right)$.
Let $\mathcal{C}=\bigcup_{A\in\mathcal{N}}\mathcal{C}_{A}$ and $\mathcal{F}$
is the smallest sigma field containing $\mathcal{C}$. Write
$\mathcal{T}_{\Lambda}=\mathcal{F}_{\mathrm{V}\backslash\Lambda}$ and
$\mathcal{T}$ for the tail- $\sigma$-algebra, i.e., intersection of
$\mathcal{T}_{\Lambda}$ over all finite subsets $\Lambda$ of $\mathbb{L}$ the
sets in $\mathcal{T}$ are called tail-measurable sets.

######  Definition 1.

[5] Let
$P_{\Lambda}:\Omega\rightarrow\overline{\mathbb{R}}:=\mathbb{R}\cup\\{-\infty,\infty\\}$
be $\mathcal{F}_{\Lambda}$-measurable mapping for all $\Lambda\in\mathcal{N}$,
then the collection $P=\left\\{P_{\Lambda}\right\\}_{\Lambda\in\mathcal{N}}$
is called a potential. Also, the following expression

$H_{\Delta,P}(\sigma)\stackrel{{\scriptstyle\text{ def
}}}{{=}}\sum_{\Delta\cap\Lambda\neq\varnothing,\Lambda\in\mathcal{N}}P_{\Lambda}(\sigma),\quad\forall\sigma\in\Omega.$
(1.1)

is called Hamiltonian $H$ associated to the potential $P$.

For a fixed inverse temperature $\beta>0$, the Gibbs specification is
determined by a family of probability kernels
$\zeta=\left(\zeta_{\Lambda}\right)_{\Lambda\in\mathcal{N}}$ defined on
$\Omega_{\Lambda}\times\mathcal{F}_{\Lambda^{c}}$ by the Boltzmann-Gibbs
weights

$\zeta_{\Lambda}(\sigma_{\Lambda}\mid\omega)=\frac{1}{\mathbf{Z}_{\Lambda}^{\omega}}e^{-\beta
H_{\Lambda,P}^{\omega}\left(\sigma_{\Lambda}\omega\right)}$ (1.2)

where $\mathbf{Z}_{\Lambda}^{\omega}=\sum_{\sigma\in\Omega_{\Lambda}}e^{-\beta
H_{\Lambda,P}^{\omega}\left(\sigma_{\Lambda}\right)}$ is the partition
function, related to free energy.

From [5], the family of mappings
$\\{\zeta_{\Lambda}(\sigma\mid\omega)\\}_{\Lambda\in\mathcal{N}}$ is the
family of proper $\mathcal{F}_{\Lambda}$-measurable quasi-probability kernels.
Thus, the collection
$\mathcal{V}=\left\\{\zeta_{\Lambda}\right\\}_{\Lambda\in\mathcal{N}}$ will be
called an $\mathbb{F}$-_specification_ if
$\zeta_{\Delta}=\zeta_{\Delta}\zeta_{\Lambda}$ whenever
$\Lambda,\Delta\in\mathcal{N}$ with $\Lambda\subseteq\Delta$. Let
$\mathcal{V}=\left\\{\zeta_{\Lambda}\right\\}_{\Lambda\in\mathcal{N}}$ be an
$\mathbb{F}$-specification; then a probability measure
$\mu\in\mathrm{P}(\mathcal{F})$ is called a _Gibbs measure with specification
$\mathcal{V}$_ if $\mu=\mu\zeta_{\Lambda}$ for each $\Lambda\in\mathcal{N}$.

##  2\. Gradient Gibbs measure

For any configuration $\omega=(\omega(x))_{x\in V}\in\mathbb{Z}^{V}$ and edge
$e=\langle x,y\rangle$ of $\vec{L}$ (oriented) the difference along the edge
$e$ is given by $\nabla\omega_{e}=\omega_{y}-\omega_{x}$ and $\nabla\omega$ is
called the gradient field of $\omega$. The gradient spin variables are now
defined by $\eta_{\langle x,y\rangle}=\omega_{y}-\omega_{x}$ for each $\langle
x,y\rangle$. The space of gradient configurations is denoted by
$\Omega^{\nabla}$. The measurable structure on the space $\Omega^{\nabla}$ is
given by $\sigma$-algebra

$\mathcal{F}^{\nabla}:=\sigma\left(\\{\eta_{e}\mid e\in\vec{L}\\}\right).$

Note that $\mathcal{F}^{\nabla}$ is the subset of $\mathcal{F}$ containing
those sets that are invariant under translations $\omega\rightarrow\omega+c$
for $c\in E$. Similarly, we define

$\mathcal{T}_{\Lambda}^{\nabla}=\mathcal{T}_{\Lambda}\cap\mathcal{F}^{\nabla},\
\mathcal{F}_{\Lambda}^{\nabla}=\mathcal{F}_{\Lambda}\cap\mathcal{F}^{\nabla}$

For nearest-neighboring (n.n.) interaction potential
$\Phi=\left(\Phi_{b}\right)_{b}$, where $b=\langle x,y\rangle$ is an edge,
define symmetric transfer matrices $Q_{b}$ by

$Q_{b}\left(\omega_{b}\right)=e^{-\left(\Phi_{b}\left(\omega_{b}\right)+|\partial
x|^{-1}\Phi_{\\{x\\}}\left(\omega_{x}\right)+|\partial
y|^{-1}\Phi_{\\{w\\}}\left(\omega_{\nu}\right)\right)}.$

Define the Markov (Gibbsian) specification as

$\gamma_{\Lambda}^{\Phi}\left(\sigma_{\Lambda}=\omega_{\Lambda}\mid\omega\right)=\left(Z_{\Lambda}^{\Phi}\right)(\omega)^{-1}\prod_{b\cap\Lambda\neq
0}Q_{b}\left(\omega_{b}\right).$

If for any bond $b=\langle x,y\rangle$ the transfer operator
$Q_{b}\left(\omega_{b}\right)$ is a function of gradient spin variable
$\zeta_{b}=\omega_{y}-\omega_{x}$ then the underlying potential $\Phi$ is
called a gradient interaction potential. Note that for all
$A\in\mathcal{F}^{\nabla}$, the kernels $\gamma_{\Lambda}^{\Phi}(A,\omega)$
are $\mathcal{F}^{\nabla}$-measurable functions of $\omega$, it follows that
the kernel sends a given measure $\mu$ on
$\left(\Omega,\mathcal{F}^{\nabla}\right)$ to another measure
$\mu\gamma_{\Lambda}^{\Phi}$ on $\left(\Omega,\mathcal{F}^{\nabla}\right)$. A
measure $\mu$ on $\left(\Omega,\mathcal{F}^{\nabla}\right)$ is called a
_gradient Gibbs measure_ if it satisfies the equality
$\mu\gamma_{\Lambda}^{\Phi}=\mu$ (detail in [10, 12, 14]).

Note that, if $\mu$ is a Gibbs measure on $(\Omega,\mathcal{F})$, then its
restriction to $\mathcal{F}^{\nabla}$ is a gradient Gibbs measure.

A boundary law is called $q$-periodic if $l_{xy}(i+q)=l_{xy}(i)$ for every
oriented edge $\langle x,y\rangle\in$ $\vec{L}$ and each $i\in\mathbb{Z}$.

It is known that there is a one-to-one correspondence between boundary laws
and tree indexed Markov chains if the boundary laws are normalisable in the
sense of Zachary [15]:

######  Definition 2.

(Normalisable boundary laws). A boundary law $l$ is said to be normalisable if
and only if

$\sum_{\omega_{x}\in\mathbb{Z}}\left(\prod_{z\in\partial
x}\sum_{\omega_{z}\in\mathbb{Z}}Q_{zx}\left(\omega_{x},\omega_{z}\right)l_{zx}\left(\omega_{z}\right)\right)<\infty$

for any $x\in V$.

The correspondence now reads the following:

######  Theorem 1.

(Theorem 3.2 in [16]). For any Markov specification $\gamma$ with associated
family of transfer matrices $\left(Q_{b}\right)_{b\in L}$ we have

1\. Each normalisable boundary law $\left(l_{xy}\right)_{x,y}$ for
$\left(Q_{b}\right)_{b\in L}$ defines a unique tree-indexed Markov chain
$\mu\in\mathcal{G}(\gamma)$ via the equation given for any connected set
$\Lambda\in\mathcal{S}$

$\mu\left(\sigma_{\Lambda\cup\partial\Lambda}=\omega_{\Lambda\cup\partial\Lambda}\right)=\left(Z_{\Lambda}\right)^{-1}\prod_{y\in\partial\Lambda}l_{yy_{\Lambda}}\left(\omega_{y}\right)\prod_{b\cap\Lambda\neq\emptyset}Q_{b}\left(\omega_{b}\right),$
(2.1)

where for any $y\in\partial\Lambda,y_{\Lambda}$ denotes the unique n.n. of $y$
in $\Lambda$.

2\. Conversely, every tree-indexed Markov chain $\mu\in\mathcal{G}(\gamma)$
admits a representation of the form (3.15) in terms of a normalisable boundary
law (unique up to a constant positive factor).

The Markov chain $\mu$ defined in (2.1) has the transition probabilities

$P_{xy}(i,j)=\mu\left(\sigma_{y}=j\mid\sigma_{x}=i\right)=\frac{l_{yx}(j)Q_{yx}(j,i)}{\sum_{s}l_{yx}(s)Q_{yx}(s,i)}.$
(2.2)

The expressions (2.2) may exist even in situations where the underlying
boundary law $\left(l_{xy}\right)_{x,y}$ is not normalisable. However, the
Markov chain given by (2.2), in general, does not have an invariant
probability measure. Therefore in $[8],[9],[10],[11]$ some nonnormalisable
boundary laws are used to give gradient Gibbs measures.

Now we give some results of above-mentioned papers. Consider a model on Cayley
tree $\Gamma^{k}=(V,\vec{L})$, where the spin takes values in the set of all
integer numbers $\mathbb{Z}$. The set of all configurations is
$\Omega:=\mathbb{Z}^{V}$.

Now we consider the following Hamiltonian:

$H(\sigma)=-J\sum_{\langle
x,y\rangle}\alpha(|\sigma_{x}-\sigma_{y}|)|\sigma_{x}-\sigma_{y}|,$ (2.3)

where

$\alpha(|m|)=\begin{cases}p,&\text{if $m\in 2\mathbb{Z}$}\\\ q,&\text{if $m\in
2\mathbb{Z}+1$}.\end{cases}$

Note that if $p=q$ then the considered model is called SOS model.

For the Hamiltonian (2.2) the transfer operator is defined by

$Q(i,j)=e^{-\beta J\alpha(|i-j|)|i-j|},$

where $\beta>0$ is the inverse temperature and $J\in\mathbb{R}$.

Also, the boundary law equation of the Hamiltonian can be written as:

$z_{i}=\left(\frac{Q(i,0)+\sum_{j\in Z_{0}}Q(i,j)z_{j}}{Q(0,0)+\sum_{j\in
Z_{0}}Q(0,j)z_{j}}\right)^{k}.$ (2.4)

Put $\theta:=\exp(-J\beta)<1$. For translation invariant boundary law, the
transfer operator $Q$ reads $Q(i-j)=\theta^{|i-j|}$ for any
$i,j\in\mathbb{Z}$. If $\theta:=e^{-J\beta}<1$ then we can write the equation
(2.4) as

$z_{i}=\left(\frac{\theta^{\alpha(|i|)|i|}+\sum_{j\in
Z_{0}}\theta^{\alpha(|i-j|)|i-j|}z_{j}}{1+\sum_{j\in
Z_{0}}\theta^{\alpha(|j|)|j|}z_{j}}\right)^{k},\
i\in\mathbb{Z}_{0}:=\mathbb{Z}\setminus\\{0\\}.$ (2.5)

Let $\\{z_{i}\\}_{i\in\mathbb{Z}}$ be $q$-periodic sequence, i.e.
$z_{i}=z_{i+q}$ for all $i\in\mathbb{Z}$.

$\left\\{\begin{array}[]{ll}z_{1}=\left(\frac{\theta^{\alpha(|1|)}+\sum_{j\in
Z_{0}}\theta^{\alpha(|1-j|)|1-j|}z_{j}}{1+\sum_{j\in
Z_{0}}\theta^{\alpha(|j|)|j|}z_{j}}\right)^{k};\\\
z_{2}=\left(\frac{\theta^{2\alpha(|2|)}+\sum_{j\in
Z_{0}}\theta^{\alpha(|2-j|)|2-j|}z_{j}}{1+\sum_{j\in
Z_{0}}\theta^{\alpha(|j|)|j|}z_{j}}\right)^{k};\\\
...\qquad...\qquad...\qquad...\qquad...\\\
z_{q}=\left(\frac{\theta^{q\alpha(|q|)}+\sum_{j\in
Z_{0}}\theta^{\alpha(|q-j|)|q-j|}z_{j}}{1+\sum_{j\in
Z_{0}}\theta^{\alpha(|j|)|j|}z_{j}}\right)^{k}.\end{array}\right.$ (2.6)

Let $u_{i}=u_{0}\sqrt[k]{z_{i}}$ for some $u_{0}>0$ then since (2.5) we obtain

$u_{i}=\frac{...+\theta^{3q}u_{i-3}^{k}+\theta^{2p}u_{i-2}^{k}+\theta^{q}u_{i-1}^{k}+u_{i}^{k}+\theta^{q}u_{i+1}^{k}+\theta^{2p}u_{i+2}^{k}+\theta^{3q}u_{i+3}^{k}+...}{...+\theta^{3q}u_{-3}^{k}+\theta^{2p}u_{-2}^{k}+\theta^{q}u_{-1}^{k}+u_{0}^{k}+\theta^{q}u_{1}^{k}+\theta^{2p}u_{2}^{k}+\theta^{3q}u_{3}^{k}+...}.$

We rewrite the last system of equations in the following form:

$u_{i}=\frac{\sum_{j=1}^{\infty}\theta^{2pj}u_{i-2j}^{k}+\sum_{j=1}^{\infty}\theta^{(2j-1)q}u_{i-2j+1}^{k}+u_{i}^{k}+\sum_{j=1}^{\infty}\theta^{(2j-1)q}u_{i+2j-1}^{k}+\sum_{j=1}^{\infty}\theta^{2pj}u_{i+2j}^{k}}{\sum_{j=1}^{\infty}\theta^{2pj}u_{-2j}^{k}+\sum_{j=1}^{\infty}\theta^{(2j-1)q}u_{-2j+1}^{k}+u_{0}^{k}+\sum_{j=1}^{\infty}\theta^{(2j-1)q}u_{2j-1}^{k}+\sum_{j=1}^{\infty}\theta^{2pj}u_{2j}^{k}},$
(2.7)

where $i\in\mathbb{Z}.$

######  Proposition 1.

Let $\\{z_{i}\\}_{i\in\mathbb{Z}}$ be $q$-periodic sequence. Then finding
$q$-periodic solutions to the system (2.5) is equivalent to solving the system
of equations (2.6).

###### Proof.

To prove the Proposition, it is sufficient to show $z_{i}=z_{q+i}$ for all
$i\in\\{1,2,...,q-1,q\\}$. Since $z_{0}=0$, for a fixed $i_{0}\in\mathbb{Z}$,
the numerator of the fraction in (2.5) can be written as

$\theta^{\alpha(|i_{0}|)|i_{0}|}+\sum_{j\in\mathbb{Z}_{0}}\theta^{\alpha(|i_{0}-j|)|i_{0}-j|}z_{j}=\sum_{j\in\mathbb{Z}}\theta^{\alpha(|i_{0}-j|)|i_{0}-j|}z_{j}.$

Also, it can be rewritten as

$\sum_{j\in\mathbb{Z}}\theta^{\alpha(|i_{0}-j|)|i_{0}-j|}z_{j}=...+\theta^{2p}z_{i_{0}-2}+\theta^{q}z_{i_{0}-1}+z_{i_{0}}+\theta^{q}z_{i_{0}+1}+\theta^{2p}z_{i_{0}+2}+...$
(2.8)

Similarly, for $i_{0}+q$ we have

$\sum_{j\in\mathbb{Z}}\theta^{\alpha(|i_{0}+q-j|)|i_{0}+q-j|}z_{j}=...+\theta^{2p}z_{i_{0}+q-2}+\theta^{q}z_{i_{0}+q-1}+z_{i_{0}+q}+\theta^{q}z_{i_{0}+q+1}+\theta^{2p}z_{i_{0}+q+2}+...$
(2.9)

If we change $z_{k+q}$ in (2.9) to $z_{k}$ for all $k\in\mathbb{N}$ then we
obtain (2.8). Namely, we have proved

$z_{i_{0}}=\left(\frac{\theta^{i_{0}\alpha(|i_{0}|)}+\sum_{j\in
Z_{0}}\theta^{\alpha(|i_{0}-j|)|i_{0}-j|}z_{j}}{i_{0}+\sum_{j\in
Z_{0}}\theta^{\alpha(|j|)|j|}z_{j}}\right)^{k}=\qquad\qquad\qquad\ \ \ \ \ \ \
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ $ $\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\ \ \ \qquad=\left(\frac{(\theta^{i_{0}+q)\alpha(|i_{0}+q|)}+\sum_{j\in
Z_{0}}\theta^{\alpha(|i_{0}+q-j|)|i_{0}+q-j|}z_{j}}{i_{0}+q+\sum_{j\in
Z_{0}}\theta^{\alpha(|j|)|j|}z_{j}}\right)^{k}=z_{i_{0}+q}.$

∎

##  3\. 4-periodic boundary laws for $k\geq 2$

In this section, we find periodic solutions (defined in [11]) to (2.5) which
correspond to periodic boundary condition. Namely, for all $m\in\mathbb{Z}$ we
consider the following sequence:

$u_{n}=\begin{cases}1,\ \emph{if }n=2m;\\\ a,\ \emph{if }n=4m-1;\\\ b,\
\emph{if }n=4m+1,\\\ \end{cases}$ (3.1)

where $a$ and $b$ are some positive numbers.

By Proposition 1, finding solutions that are formed in (3.1) to (2.5) is
equivalent to solving the following system of equations:

$a=\frac{...+\theta^{4p}a^{k}+\theta^{3q}+\theta^{2p}b^{k}+\theta^{q}+a^{k}+\theta^{q}+\theta^{2p}b^{k}+\theta^{3q}+\theta^{4p}a^{k}+...}{...+\theta^{4p}+\theta^{3q}b^{k}+\theta^{2p}+\theta^{q}a^{k}+1+\theta^{q}b^{k}+\theta^{2p}+\theta^{3q}a^{k}+\theta^{4p}+...};$

$b=\frac{...+\theta^{4p}b^{k}+\theta^{3q}+\theta^{2p}a^{k}+\theta^{q}+b^{k}+\theta^{q}+\theta^{2p}a^{k}+\theta^{3q}+\theta^{4p}b^{k}+...}{...+\theta^{4p}+\theta^{3q}b^{k}+\theta^{2p}+\theta^{q}a^{k}+1+\theta^{q}b^{k}+\theta^{2p}+\theta^{3q}a^{k}+\theta^{4p}+...}.$

Namely,

$a=\frac{2(\theta^{q}+\theta^{3q}+...)+(1+2\theta^{4p}+2\theta^{8p}+...)a^{k}+2(\theta^{2p}+\theta^{6p}+...)b^{k}}{1+2\theta^{2p}+2\theta^{4p}+...+(\theta^{q}+\theta^{3q}+...)(a^{k}+b^{k})};$

$b=\frac{2(\theta^{q}+\theta^{3q}+...)+(1+2\theta^{4p}+2\theta^{8p}+...)b^{k}+2(\theta^{2p}+\theta^{6p}+...)a^{k}}{1+2\theta^{2p}+2\theta^{4p}+...+(\theta^{q}+\theta^{3q}+...)(a^{k}+b^{k})}.$

Taking into account $\theta<1$ one writes the last system of equations as
follows:

$a=\frac{\frac{2\theta^{q}}{1-\theta^{2q}}+\frac{1+\theta^{4p}}{1-\theta^{4p}}a^{k}+\frac{2\theta^{2p}}{1-\theta^{4p}}b^{k}}{\frac{1+\theta^{2p}}{1-\theta^{2p}}+\frac{\theta^{q}}{1-\theta^{2q}}(a^{k}+b^{k})},\quad
b=\frac{\frac{2\theta^{q}}{1-\theta^{2q}}+\frac{1+\theta^{4p}}{1-\theta^{4p}}b^{k}+\frac{2\theta^{2p}}{1-\theta^{4p}}a^{k}}{\frac{1+\theta^{2p}}{1-\theta^{2p}}+\frac{\theta^{q}}{1-\theta^{2q}}(a^{k}+b^{k})}.$
(3.2)

For all $k\in\mathbb{N}$ and $p,q\in\mathbb{R}$, the analysis of (3.2) is so
difficult and that’s why we consider the case $2p=q=1$ and $k=2.$ Then (3.2)
can be written as

$a=\frac{\tau a^{2}+2b^{2}+2}{a^{2}+b^{2}+\tau+2},\quad b=\frac{\tau
b^{2}+2a^{2}+2}{a^{2}+b^{2}+\tau+2},$ (3.3)

where $\tau=\theta+\frac{1}{\theta}>2.$

At first, we consider the case $a=b$. The system of equations (3.3) is reduced
to the polynomial equation:

$2a^{3}-(\tau+2)a^{2}+(\tau+2)a-2=0.$ (3.4)

Since the last equation has a solution $a=1$, we divide both sides of (3.4) by
$a-1$. Consequently, one gets

$2a^{2}-\tau a+2=0.$

If $\tau>4$ then the last quadratic equation has two solutions.

######  Theorem 2.

Let $\tau=J\beta+\frac{1}{J\beta}$. Then for the model (2.3) on the the Cayley
tree of order two the following assertions hold:

  1. (1)

If $\tau\leq 4$, then there is precisely one GGM associated to a 2-periodic
boundary law.

  2. (2)

If $\tau>4$, then there are precisely three such GGMs.

It is important to consider Gradient Gibbs measures associated with a
4-periodic boundary law for the model (2.3). The following theorem gives us a
full description of Gradient Gibbs measures associated to a 4-periodic
boundary law.

######  Theorem 3.

Let $\tau=J\beta+\frac{1}{J\beta}$, $\tau_{cr}^{(1)}\approx 5.73$, and
$\tau_{cr}^{(2)}\approx 6.261$. Then for Gradient Gibbs measures associated
with a 4-periodic boundary law for the model (2.3) on the Cayley tree of order
two the following statements hold:

  1. (1)

If $\tau\leq\tau_{cr}^{(1)}$, then there is not any GGM.

  2. (2)

If $\tau_{cr}^{(1)}<\tau<\tau_{cr}^{(2)}$, then there exists a unique GGM.

  3. (3)

If $\tau=\tau_{cr}^{(2)}$, then there are exactly two GGMs.

  4. (4)

If $\tau_{cr}^{(2)}<\tau<8$, then there are exactly three GGMs.

  5. (5)

If $\tau=8$, then there are exactly four GGMs.

  6. (6)

If $\tau>8$, then there are exactly five such GGMs.

###### Proof.

Now we consider the case $a\neq b$. Then the system of equations (3.3) can be
written as

$\begin{cases}a^{3}+ab^{2}+(\tau+2)a=\tau a^{2}+2b^{2}+2;\\\\[5.69054pt]
b^{3}+a^{2}b+(\tau+2)b=\tau b^{2}+2a^{2}+2.\end{cases}$ (3.5)

Now we subtract the second equation of (3.5) from the first one and get

$a^{3}-b^{3}-ab(a-b)+(\tau+2)(a-b)=(\tau-2)(a^{2}-b^{2}).$

Since $a\neq b$, both sides can be divided by $a-b$ and one gets

$a^{2}+b^{2}+\tau+2=(\tau-2)(a+b).$ (3.6)

By adding the second and first equations of (3.5), we have

$(\tau-2)ab=2(\tau-2)(a+b)-2(\tau+1)$ (3.7)

Let $a+b=x$ and $ab=y$. By using (3.6) and (3.7) one gets a new system of
equations with respect to $x$ and $y$ that is equivalent to (3.5):

$\begin{cases}x^{2}-2y+\tau+2=(\tau-2)x\\\
(\tau-2)y=2(\tau-2)x-2(\tau+1)\end{cases}$ (3.8)

In order to find the number of solutions of the last system we can consider
the following quadratic equation with respect to $x$:

$(\tau-2)x^{2}-(\tau^{2}-4)x+\tau^{2}+4\tau=0.$ (3.9)

It is easy to check that

$x_{1}=\frac{(\tau-2)(\tau+2)-\sqrt{(\tau-2)(\tau^{3}-2\tau^{2}-20\tau-8)}}{2(\tau-2)}$

and

$x_{2}=\frac{(\tau-2)(\tau+2)+\sqrt{(\tau-2)(\tau^{3}-2\tau^{2}-20\tau-8)}}{2(\tau-2)}$

are solutions to the equation (3.9). Put
$P(\tau)=\tau^{3}-2\tau^{2}-20\tau-8$. From $\tau>0$ it is sufficient to find
only positive roots of $P(\tau)$. By Descartes’ theorem (e.g. [13]) $P(\tau)$
has at most one positive root. On the other hand, $P(5)<0$ and $P(6)>0$ i.e.,
by Intermediate Value Theorem $P(\tau)$ has at least one root in the segment
$[5,6]$. Hence $P(\tau)$ has exactly one positive root. Let
$\tau_{cr}(\tau_{cr}^{(1)}\approx 5.73)$ be the positive root of the
polynomial. Consequently, we can conclude that if $2<\tau<\tau_{cr}^{(1)}$
then the system of equations (3.8) has not any positive solution. Let
$\tau=\tau^{(1)}_{cr}$, then the system (3.8) has exactly one positive root.
For the case $\tau>\tau_{cr}^{(1)}$, then (3.8) has exactly two positive roots
if we can show $x_{1}>0$. Namely, after short calculations, if
$\tau>\tau_{cr}^{(1)}$ then we can show the inequality

$\frac{(\tau-2)(\tau+2)-\sqrt{(\tau-2)(\tau^{3}-2\tau^{2}-20\tau-8)}}{2(\tau-2)}>0$

is equivalent to the inequality $\tau^{2}+4\tau>0$.

For the case $\tau>\tau_{cr}^{(1)}$, from $a+b=x_{i}$ and $ab=y_{i}$
($i\in\\{1,2\\}$) after short calculations, we have two quadratic equations
respectively to $x_{1}$ and $x_{2}$:

$2(\tau-2)a^{2}-(\tau^{2}-4-\sqrt{(\tau-2)(\tau^{3}-2\tau^{2}-20\tau-8)})a+\ \
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\ \ \ $ $\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\ \ \ \ \ \ \ \ \
+2(\tau^{2}-4-\sqrt{(\tau-2)(\tau^{3}-2\tau^{2}-20\tau-8)}-2(\tau+1))=0$

and

$2(\tau-2)a^{2}-(\tau^{2}-4+\sqrt{(\tau-2)(\tau^{3}-2\tau^{2}-20\tau-8)})a+\ \
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\ \ \ $ $\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\ \ \ \ \ \ \ \ \
+2(\tau^{2}-4+\sqrt{(\tau-2)(\tau^{3}-2\tau^{2}-20\tau-8)}-2(\tau+1))=0.$

The discriminants are

$D_{1,2}(\tau)=2(\tau-2)(\tau^{3}-8\tau^{2}+4\tau+40\pm(6-\tau)\sqrt{(\tau-2)(\tau^{3}-2\tau^{2}-20\tau-8)}).$

Now we find positive zeroes of

$\tau^{3}-8\tau^{2}+4\tau+40\pm(6-\tau)\sqrt{(\tau-2)(\tau^{3}-2\tau^{2}-20\tau-8)}=0.$

A solution to the last equation is also the solution to the following
equation:

$\tau^{4}-16\tau^{3}+80\tau^{2}-160\tau+256=(\tau-8)\left(\tau^{3}-8\tau^{2}+16\tau-32\right)=0.$

The discriminant of the polynomial $Q(\tau):=\tau^{3}-8\tau^{2}+16\tau-32$ is
negative and that is why $Q(\tau)$ has a unique real root. Since $Q(0)<0$ and
$\lim_{\tau\rightarrow\infty}Q(\tau)=+\infty$ this root $\tau_{cr}^{(2)}$ is
positive and $\tau_{cr}^{(2)}\approx 6.261$. Consequently, 8 (resp.
$\tau_{cr}^{(2)}$) is a unique real positive root of $D_{1}(\tau)$
($D_{2}(\tau)$).

Finally, we consider the case $\tau=\tau^{(1)}_{cr}$. From above, it is
sufficient to solve the following equation:

$2(\tau-2)a^{2}-(\tau^{2}-4)a+2(\tau^{2}-2\tau-6))=0.$

Its discriminant is

$D(\tau)=2(\tau-2)(\tau^{3}-8\tau^{2}+4\tau+40).$

It’s easy to check $D(\tau^{(1)}_{cr})<0$, thus there no positive solution to
(3.8). ∎

## Acknowledgements

The work supported by the fundamental project (number: F-FA-2021-425) of The
Ministry of Innovative Development of the Republic of Uzbekistan. I thank
Professor U.A.Rozikov for useful discussions and suggestions which have
improved the paper.

## Statements and Declarations

Conflict of interest statement: The author states that there is no conflict of
interest.

## Data availability statements

The datasets generated during and/or analyzed during the current study are
available from the corresponding author upon reasonable request.

## References

  * [1] M. Biskup, R. Kotecky, Phase coexistence of gradient Gibbs states, Probab. Theory Related Fields 139(1-2) (2007), 1-39.
  * [2] R. Bissacot, E.O. Endo, A.C.D. van Enter, Stability of the phase transition of critical-field Ising model on Cayley trees under inhomogeneous external fields, Stoch. Process. Appl. 127(12) (2017), 4126-4138.
  * [3] G.I. Botirov, F.H. Haydarov, Gradient Gibbs measures for the SOS model with integer spin values on a Cayley tree. J. Stat. Mech. Theory Exp. 9 (2020), 093102, 9 pp.
  * [4] A.C.D. van Enter, C.Kulske, Non-existence of random gradient Gibbs measures in continuous interface models in $d=2$, Ann. Appl. Probab. 18 (2008) 109-119.
  * [5] S. Friedli, Y. Velenik, Statistical mechanics of lattice systems. A concrete mathematical introduction. Cambridge University Press, Cambridge, 2018. xix +622 pp.
  * [6] T. Funaki, H. Spohn, Motion by mean curvature from the Ginzburg-Landau $\nabla\varphi$ interface model, Comm. Math. Phys., 185, 1997.
  * [7] H.O. Georgii, Gibbs Measures and Phase Transitions, Second edition. de Gruyter Studies in Mathematics, 9. Walter de Gruyter, Berlin, 2011.
  * [8] F. Henning, C. Kulske, A. Le Ny, U.A. Rozikov, Gradient Gibbs measures for the SOS model with countable values on a Cayley tree. Electron. J. Probab. 24 (2019), Paper No. 104, 23 pp.
  * [9] F. Henning, C. Kulske, Existence of gradient Gibbs measures on regular trees which are not translation invariant. Ann. Appl. Probab. 33(4): 3010-3038, (2023).
  * [10] F. Henning, C. Kulske, Coexistence of localized Gibbs measures and delocalized gradient Gibbs measures on trees. Ann. Appl. Probab. 31(5): 2284-2310, (2021).
  * [11] F.H. Haydarov, U.A. Rozikov, Gradient Gibbs measures of an SOS model on Cayley trees: 4-periodic boundary laws, Reports on Mathematical Physics, 90(1), (2022), 81-101.
  * [12] C. Kulske, P. Schriever. Gradient Gibbs measures and fuzzy transformations on trees, Markov Process. Relat. Fields 23, (2017), 553-590.
  * [13] V.V. Prasolov, Polynomials. Springer, Berlin. (2004).
  * [14] U.A. Rozikov, Gibbs measures on Cayley trees. World Sci. Publ. Singapore. (2013).
  * [15] S. Sheffield, Random surfaces: Large deviations principles and gradient Gibbs measure classifications. Thesis (Ph.D.)-Stanford University, (2003), 205 pages.
  * [16] S. Zachary, Countable state space Markov random fields and Markov chains on trees, Ann. Probab. 11(4) (1983), 894-903.