Self-excited waves in complex social systems

V.I. Yukalov1,2,∗, E.P. Yukalova3

1Bogolubov Laboratory of Theoretical Physics,
Joint Institute for Nuclear Research, Dubna 141980, Russia

and

Instituto de Fisica de São Carlos, Universidade de São Paulo,
CP 369, São Carlos 13560-970, São Paulo, Brazil

3Laboratory of Information Technologies,
Joint Institute for Nuclear Research, Dubna 141980, Russia

∗Corresponding author e-mail<EMAIL_ADDRESS>

###### Abstract

A social system is considered whose agents choose between several alternatives
of possible actions. The system is described by the fractions of agents
preferring the corresponding alternatives. The agents interact with each other
by exchanging information on their choices. Each alternative is characterized
by three attributes: utility, attractiveness, and replication. The agents are
heterogeneous having different initial conditions and different types of
memory, which can be long-term or short-term. The agent interactions,
generally, can depend on the distance between the agents, varying from short-
range to long-range interactions. The emphasis in the paper is on long-range
interactions. In a mixed society consisting of agents with both long-term and
short-term memory, there appears the effect of spontaneous excitation of
preference waves, when the fractions of agents preferring this or that
alternative suddenly start strongly oscillating, either periodically or
chaotically. Since the considered society forms a closed system, without any
external influence, the arising waves are internally self-excited by the
society.

Keywords: complex social systems, information exchange, long-term memory,
short-term memory, herding effect, spontaneous preference waves

##  1 Introduction

Understanding the behaviour of social systems is of long standing and
permanent interest [1], and it has been described employing various models
based on statistical physics (see, e.g., reviews [2, 3]) or on the theory of
networks [5, 6, 7, 8, 9, 10]. The recent review [4] summarizes both
approaches. Here we keep in mind social networks consisting of agents choosing
between different alternatives. Interactions between society members usually
are described by given functions, similarly to the interactions between
particles or spins in statistical mechanics. In that sense, society members
are treated as nodes mechanically accomplishing prescribed actions. In the
case of realistic biological, especially human, societies, such models,
clearly, give a rather simplified picture, since the agents of these societies
are not mechanical devices, but are intelligent agents. The basic feature of
an intelligent agent is the ability, after evaluating the available
information, to make decisions choosing between several alternative actions
[11, 12, 13, 14, 15].

The aim of the present paper, is to suggest a model of a society composed of
intelligent agents. To our understanding, such a model could provide a more
realistic description of social systems consisting of intelligent agents
taking decisions with respect to available alternatives. The basic points of
the suggested theory are as follows.

(i) The approach is probabilistic. Each agent is characterized by the
probability of choosing this or that alternative. This probability has the
meaning of a frequentist measure showing the fraction of agents choosing the
given alternative.

(ii) The probability measure takes into account three factors: the utility of
alternatives, their attractiveness, and the attitude of agents towards
replicating the actions of other members of the society. Thus, in addition to
the estimation of utility of alternatives, our model includes the influence of
emotions, described by the attractiveness of alternatives, and takes into
account the herding effect.

Endeavors of including emotions in the process of choice have been undertaken
in the frame of quantum decision theory [16, 17, 18, 19, 20, 21, 22, 23].
However, as has been shown [24], the use of quantum theory can be avoided, and
an approach can be developed incorporating emotions into decision process,
without resorting to quantum techniques. Below we suggest an approach of
describing the society of intelligent agents, employing only classical
notions. We study a novel kind of complex society composed of agents with
different types of memory, so that a fraction of the society members is
endowed with long-term memory, while the others have short-term memory. In
this complex society, interesting new effects appear, such as suddenly arising
self-excited oscillations of preferences. These spontaneous oscillations can
be either periodic or chaotic.

The layout of the paper is as follows. In Sec. 2, we formulate the model of a
society of intelligent agents deciding between several alternatives. The
choice is based of three attributes, utility, attractiveness, and herding. In
Sec. 3, the problem is specified to the consideration of two groups of agents
choosing between two alternatives. The groups are composed of agents
possessing different types of memory, either long-term or short-term memory.
Section 4 presents a detailed investigation of different dynamic regimes that
can occur in the process of decision making. Section 5 concludes.

##  2 Society of intelligent agents

Consider a society of $N$ agents, enumerated by the index $j=1,2,\ldots,N$.
Assume that each agent needs to choose between $N_{A}$ alternative actions
$A_{n}$, numbered by the index $n=1,2,\ldots,N_{A}$. This can be the choice
between several candidates at elections, between different goods in a shop,
between several jobs, etc.

The probability that a $j$-th agent chooses an alternative $A_{n}$ at time $t$
is denoted as $p_{j}(A_{n},t)$. By the meaning of probability,
$p_{j}(A_{n},t)$ has to be non-negative and normalized,

$\sum_{n=1}^{N_{A}}p_{j}(A_{n},t)=1\;,\qquad 0\leq p_{j}(A_{n},t)\leq 1\;.$
(1)

The utility of an alternative $A_{n}$ defines the weight $f(A_{n},t)$ ascribed
to this alternative, because of which $f(A_{n},t)$ can be called utility
factor. Being a weight, it enjoys the properties of classical probability,

$\sum_{n=1}^{N_{A}}f_{j}(A_{n},t)=1\;,\qquad 0\leq f_{j}(A_{n},t)\leq 1\;.$
(2)

The explicit expression for the utility factor can be found from the
minimization of an information functional [20, 25].

Important role in the process of decision making is played by emotions that
can be characterized by the attraction factor $q(A_{n},t)$. Emotions can be
positive or negative, because of which the attraction factor varies in the
interval $[-1,1]$,

$-1\leq q_{j}(A_{n},t)\leq 1\;.$ (3)

If the agents of a society exchange information with each other, there appears
herding effect when the members of the society incline to replicate the
actions of others. Let this effect be quantified by the herding factor
$h(A_{n},t)$, being in the interval

$-1\leq h_{j}(A_{n},t)\leq 1\;.$ (4)

The probability of preferring an alternative $A_{n}$ is the superposition of
the utility, attraction, and herding factors. Taking into account the fact
that the process of making a decision requires some time, say $\tau$, we have

$p_{j}(A_{n},t+\tau)=f_{j}(A_{n},t)+q_{j}(A_{n},t)+h_{j}(A_{n},t)\;.$ (5)

This expression differentiates the probabilistic approach we follow from
multi-attribute utility theory, where the expected utility functional is
represented as a superposition of weighted parts associated with different
attributes [26, 27].

From the normalization conditions (1) and (2), it follows that the sum of the
attraction factor and herding factor over all alternatives is zero:

$\sum_{n}\left[\;q_{j}(A_{n},t)+h_{j}(A_{n},t)\;\right]=0\;.$ (6)

Assuming that attraction and herding are independent, we have

$\sum_{n}q_{j}(A_{n},t)=0\;,\qquad\sum_{n}h_{j}(A_{n},t)=0\;.$ (7)

The latter equations can be called alternation law.

Although, on a long time scale there exists time discounting [28], however the
utility of alternatives varies in time slowly, which allows us to consider it
constant,

$f_{j}(A_{n},t)=f_{j}(A_{n})\;.$ (8)

On the contrary, information exchange between the society agents is a fast
process, such that the attraction factor of a $j$-th agent essentially depends
on the amount of information $M_{j}(t)$ obtained by this agent. The attraction
factor can be modeled [29] by the expression

$q_{j}(A_{n},t)=q_{j}(A_{n})\exp\\{-M_{j}(t)\\}\;,$ (9)

where $q(A_{n})=q(A_{n},0)$ is an initial value of the attraction factor. This
expression is similar to the time dependence of decoherence factor under
nondestructive repeated measurements [30, 31]. The amount of information
obtained by a $j$-th agent by the time $t$ can be written as

$M_{j}(t)=\sum_{t^{\prime}=0}^{t}\;\sum_{i=1}^{N}J_{ji}(t,t^{\prime})\;\mu_{ji}(t^{\prime})\;.$
(10)

This can be called remembered information. Here $J_{ji}(t,t^{\prime})$ is the
intensity of interactions during the information transfer from an $i$-th agent
to the $j$-th agent in the period of time between $t^{\prime}$ and $t$. The
information gain, received by the $j$-th agent from an $i$-th agent at time
$t$, is given in the form of the Kullback-Leibler [32, 33] relative
information

$\mu_{ji}(t)=\sum_{n=1}^{N_{A}}p_{j}(A_{n},t)\;\ln\;\frac{p_{j}(A_{n},t)}{p_{i}(A_{n},t)}\;.$
(11)

Note that the information gain enjoys the properties $\mu_{ji}\geq 0$ and
$\mu_{jj}=0$. At the initial moment of time, no information has yet been
transferred, implying that

$M_{j}(0)=0\;.$ (12)

The interactions $J_{ij}(t,t^{\prime})$, in general, depend on the distance
between the agents. In the case of short-range interactions the topology of
the agent locations is important. In the opposite case of long-range
interactions, the geometry of the agent network plays no role, with the
interactions acquiring the form

$J_{ij}(t,t^{\prime})=\frac{1}{N-1}\;J(t,t^{\prime})\;.$ (13)

Our concern here is a society composed of intelligent agents, like humans who
are able to interact with each other irrespectively of the distance between
them. This kind of distance-independent interactions are provided nowadays by
internet, mass media, and phones. In addition, the agents are not located at
fixed nodes but can freely move varying their whereabouts. Keeping in mind
these location-independent interactions, we accept in what follows the form
(13). Then the amount of information kept in the agent memory reads as

$M_{j}(t)=\sum_{t^{\prime}=0}^{t}\frac{J(t,t^{\prime})}{N-1}\;\sum_{i=1}^{N}\mu_{ji}(t^{\prime})\;.$
(14)

The agent memory is characterized by its duration. In two opposite situations,
it can be long-term or short-term. The ultimate long-term memory is permanent
in time,

$J(t,t^{\prime})=J\qquad(long-term)\;.$ (15)

Then the amount of remembered information (14) takes the form

$M_{j}(t)=\frac{J}{N-1}\;\sum_{t^{\prime}=0}^{t}\;\sum_{i=1}^{N}\mu_{ji}(t^{\prime})\;.$
(16)

The opposite case is the short-term memory, when the information from only the
last step is remembered,

$J(t,t^{\prime})=J\;\delta_{tt^{\prime}}\qquad(short-term)\;.$ (17)

In that case, the available remembered information (14) becomes

$M_{j}(t)=\frac{J}{N-1}\sum_{i=1}^{N}\mu_{ji}(t)\;.$ (18)

In complex societies, there exists a collective effect, when the agents are
prone to imitate the actions of others. This is called the herd effect. Herd
behavior occurs in animal herds, packs, bird flocks, fish schools and so on,
as well as in humans. It is well known and studied for many years [34, 35, 36,
37, 38, 39, 40]. In evolution equations of social and biological systems, the
mathematical description of herding is represented by the replication term
[41], which for our case takes the form

$h_{j}(A_{n},t)=\varepsilon_{j}\;\left\\{\frac{1}{N-1}\sum_{i(\neq
j)}^{N}[\;f_{i}(A_{n})+q_{i}(A_{n},t)\;]-[\;f_{j}(A_{n})+q_{j}(A_{n},t)\;]\right\\}\;.$
(19)

The parameters $\varepsilon_{j}$ describe the intensity of the herding
behaviour. Due to the normalization conditions (1), (2), (3), (4), and (7),
these parameters satisfy the inequalities

$0\leq\varepsilon_{j}\leq 1\qquad(j=1,2,\ldots,N)\;.$ (20)

Without the loss of generality, time can be measured in units of $\tau$.
Substituting the herding term (19) into expression (5) yields the equation

$p_{j}(A_{n},t+1)=(1-\varepsilon_{j})\;[\;f_{j}(A_{n})+q_{j}(A_{n},t)\;]+\frac{\varepsilon_{j}}{N-1}\sum_{i(\neq
j)}^{N}[\;f_{i}(A_{n})+q_{i}(A_{n},t)\;]\;.$ (21)

This it the evolution equation for the probability that the $j$-th agent
chooses at time $t+1$ the alternative $A_{n}$. The right-hand side of this
equation at time $t=0$ gives the initial condition

$p_{j}(A_{n},0)=(1-\varepsilon_{j})\;[\;f_{j}(A_{n})+q_{j}(A_{n},0)\;]+\frac{\varepsilon_{j}}{N-1}\sum_{i(\neq
j)}^{N}[\;f_{i}(A_{n})+q_{i}(A_{n},0)\;]\;.$ (22)

##  3 Agents with different types of memory

Now we need to specify what kind of agents we consider. Agents can differ by
initial conditions or by the types of their memory. The simplest case is when
all agents possess the same type of memory, although even then the dynamics of
preferences can exhibit rather nontrivial behaviour [29]. Here we study a more
interesting case, where the society consists of agents with different types of
memory. A fraction of agents has long-term memory and the other part, short-
term memory. Such a mixture of agents with different properties much better
models the real human societies. It turns out that this more complex society
exhibits unusual phenomena that are absent in homogeneous societies where all
agents possess the same type of memory. For instance, there appear self-
excited waves of preferences, spontaneously arising without any external
influence. These waves can appear at the beginning of dynamics or may not
exist at the beginning of decision processes, but appear suddenly after some
time of quite smooth dynamics.

Suppose the society consists of two types of agents. One part enjoys long-term
memory associated with the information function (16), while the other part has
short-term memory corresponding to the information function (18). A group of
similar agents can be represented by a frequentist probability showing the
fraction of agents choosing an alternative $A_{n}$ at time $t$. For the case
of two groups, we have the probabilities $p_{1}(A_{n},t)$ and
$p_{2}(A_{n},t)$.

Let us also consider the very often met situation when the choice is between
two alternatives, say $A_{1}$ and $A_{2}$. Then, keeping in mind the
normalization conditions (1), (2), and (7), with $j=1,2$, we can simplify the
notation for the probabilities,

$p_{j}(A_{1},t)\equiv p_{j}(t)\;,\qquad p_{j}(A_{2},t)=1-p_{j}(t)\;,$ (23)

utility factors,

$f_{j}(A_{1})\equiv f_{j}\;,\qquad f_{j}(A_{2})=1-f_{j}\;,$ (24)

attraction factors,

$q_{j}(A_{1},t)\equiv q_{j}(t)\;,\qquad q_{j}(A_{2},t)=-q_{j}(t)\;,$ (25)

and herding factors,

$h_{j}(A_{1},t)\equiv h_{j}(t)\;,\qquad h_{j}(A_{2},t)=-h_{j}(t)\;.$ (26)

Then we have the probabilities

$p_{j}(t+1)=f_{j}+q_{j}(t)+h_{j}(t)\qquad(j=1,2)\;,$ (27)

with the attraction and herding factors

$q_{j}(t)=q_{j}\exp\\{-M_{j}(t)\\}\qquad(j=1,2)\;,$
$h_{j}(t)=\varepsilon_{j}\;[\;f_{i}+q_{i}(t)-f_{j}-q_{j}(t)\;]\qquad(i\neq
j)\;,$ (28)

where $q_{j}\equiv q_{j}(0)$.

Let us mark the group of agents with long-term memory as the group number one,
and the group of agents with short-term memory, as the second group. For the
corresponding available remembered information, setting $J=1$, we have in the
case of long-term memory

$M_{1}(t)=\sum_{t^{\prime}=0}^{t}\mu_{12}(t^{\prime})$ (29)

and for short-term memory,

$M_{2}(t)=\mu_{21}(t)\;.$ (30)

The information gain reads as

$\mu_{ij}(t)=p_{i}(t)\;\ln\;\frac{p_{i}(t)}{p_{j}(t)}+[\;1-p_{i}(t)\;]\;\ln\;\frac{1-p_{i}(t)}{1-p_{j}(t)}\;,$
(31)

where $i\neq j$. Thus we come to the evolution equations

$p_{1}(t+1)=(1-\varepsilon_{1})\;[\;f_{1}+q_{1}(t)\;]+\varepsilon_{1}\;[\;f_{2}+q_{2}(t)\;]\;,$
$p_{2}(t+1)=(1-\varepsilon_{2})\;[\;f_{2}+q_{2}(t)\;]+\varepsilon_{2}\;[\;f_{1}+q_{1}(t)\;]\;,$
(32)

with the initial conditions

$p_{1}(0)=(1-\varepsilon_{1})\;[\;f_{1}+q_{1}\;]+\varepsilon_{1}\;[\;f_{2}+q_{2}\;]\;,$
$p_{2}(0)=(1-\varepsilon_{2})\;[\;f_{2}+q_{2}\;]+\varepsilon_{2}\;[\;f_{1}+q_{1}\;]\;$
(33)

that are defined by the parameters $f_{1}$, $q_{1}$, $f_{2}$, and $q_{2}$.

Recall that $p_{1}(t)$ is the fraction of agents with long-term memory
preferring the alternative $A_{1}$ at time $t$, while $p_{2}(t)$ is the
fraction of agents with short-term memory preferring the alternative $A_{1}$
at time $t$. Because of the exchange of information and the herding effect,
these fractions vary in time.

##  4 Dynamics of preferences

Before going to the numerical investigation of Eqs. (32), it is possible to
make general conclusions on the expected role of the herding effect. To this
end, it is easy to notice that, when the herding behavior is absent, hence
$\varepsilon_{1}=\varepsilon_{2}=0$, then Eqs. (32) read as

$p_{1}(t+1)=f_{1}+q_{1}(t)\qquad(\varepsilon_{1}=0)\;,$
$p_{2}(t+1)=f_{2}+q_{2}(t)\qquad(\varepsilon_{2}=0)\;.$ (34)

This does not mean that the first and second groups are independent, since
they remain connected through the information exchange entering the remembered
information (29) and (30). However, to some extent, the evolution of group
preferences can be quite different because of so different memory properties
of the groups.

Increasing the herding parameters $\varepsilon_{j}>0$ switches on the herding
behavior that, due to the replication terms $h_{j}(t)$, should smooth the
differences in the dynamics of the probabilities $p_{j}(t)$. When the herding
parameters reach the value $\varepsilon_{1}=\varepsilon_{2}=1/2$, the
behaviour of both groups becomes identical, since

$p_{1}(t+1)=p_{2}(t+1)\qquad\left(\varepsilon_{1}=\varepsilon_{2}=\frac{1}{2}\right)\;.$
(35)

With the following increase of the herding parameters, the difference in the
behaviour of the groups starts growing, reaching the maximum for
$\varepsilon_{1}=\varepsilon_{2}=1$, when the probabilities become

$p_{1}(t+1)=f_{2}+q_{2}(t)\qquad(\varepsilon_{1}=1)\;,$
$p_{2}(t+1)=f_{1}+q_{1}(t)\qquad(\varepsilon_{2}=1)\;.$ (36)

Thus the expected behaviour of the group probabilities can be rather different
at the point $\varepsilon_{1}+\varepsilon_{2}=0$, when there is no herding
effect. The difference smooths when approaching the line
$\varepsilon_{1}+\varepsilon_{2}=1$. An then the difference increases again
when drifting to the line $\varepsilon_{1}+\varepsilon_{2}=2$. Numerical
solution of Eqs. (32) demonstrates that there occur various types of
qualitatively different behaviour. In the first four Figures 1 to 4, we show
dynamic regimes corresponding to what can be called moderate herding, with the
herding parameters varying from the point $\varepsilon_{1}+\varepsilon_{2}=0$
to the line $\varepsilon_{1}+\varepsilon_{2}=1$. Figures 5 to 10 present
dynamic regimes for large herding parameters between the lines
$\varepsilon_{1}+\varepsilon_{2}=1$ and $\varepsilon_{1}+\varepsilon_{2}=2$.
Overall, the following dynamic regimes can happen.

1. In the absence of herding, the functions $p_{j}(t)$ tend to fixed points

$p_{j}^{*}=\lim_{t\rightarrow\infty}p_{j}(t)\;.$ (37)

The appearance of herding makes the trajectories closer to each other. If
there are sharp variations of the functions, they become smoothed by the
herding effect, as is shown in Fig. 1.

Figure 1: Probabilities $p_{1}(t)$ (solid line) and $p_{2}(t)$ (dash-dotted
line) for the initial conditions $f_{1}=0.4$, $f_{2}=0.1$, $q_{1}=0.59$, and
$q_{2}=0.6$. (a) $\varepsilon_{1}=\varepsilon_{2}=0$. The fixed points are
$p_{1}^{*}=0.4$ and $p_{2}^{*}=0.636$; (b) $\varepsilon_{1}=0.5$ and
$\varepsilon_{2}=0$. The consensual fixed point is $p_{1}^{*}=p_{2}^{*}=0.7$;
(c) $\varepsilon_{1}=0.5$ and $\varepsilon_{2}=0.3$. The consensual fixed
point is $p_{1}^{*}=p_{2}^{*}=0.7$.

2. Without herding, the fraction $p_{1}(t)$ of agents with long-term memory tends to a fixed point without sharp variations, while the fraction $p_{2}(t)$ of agents with short-term memory experiences at the beginning essential oscillations that attenuate and finally the function tends to a fixed point, as is demonstrated in Fig. 2.

Figure 2: Probabilities $p_{1}(t)$ (solid line) and $p_{2}(t)$ (dash-dotted
line) for the initial conditions $f_{1}=0.8$, $f_{2}=0.9$, $q_{1}=0.19$, and
$q_{2}=-0.8$. (a) $\varepsilon_{1}=\varepsilon_{2}=0$. The fixed points are
$p_{1}^{*}=0.8$ and $p_{2}^{*}=0.377$; (b) $\varepsilon_{1}=0.05$ and
$\varepsilon_{2}=0$. The fixed points are $p_{1}^{*}=0.778$ and
$p_{2}^{*}=0.362$; (c) $\varepsilon_{1}=0.1$ and $\varepsilon_{2}=0.85$. The
fixed points are $p_{1}^{*}=0.730$ and $p_{2}^{*}=0.695$.

3. When there is no herding, $p_{1}(t)$ tends to a fixed point, while $p_{2}(t)$ permanently oscillates. Switching on the herding parameters results in the permanent oscillation of both functions $p_{1}(t)$, as well as $p_{2}(t)$. Increasing the herding parameters further, after initial oscillations, forces both functions to tend to their fixed points. This is illustrated in Fig. 3.

Figure 3: Probabilities $p_{1}(t)$ (solid line) and $p_{2}(t)$ (dash-dotted
line) for the initial conditions $f_{1}=0.8$, $f_{2}=1$, $q_{1}=0.1$, and
$q_{2}=-0.99$. (a) $\varepsilon_{1}=\varepsilon_{2}=0$. The probability for
the agents with long-term memory tends to the fixed point $p_{1}^{*}=0.8$,
while the function $p_{2}(t)$ oscillates for all times; (b)
$\varepsilon_{1}=0.1$ and $\varepsilon_{2}=0$. Both functions $p_{1}(t)$ and
$p_{2}(t)$ permanently oscillate; (c) $\varepsilon_{1}=\varepsilon_{2}=0.1$.
The functions, after experiencing at the beginning oscillations, tend to the
fixed points $p_{1}^{*}=0.75$ and $p_{2}^{*}=0.353$; (d) $\varepsilon_{1}=0$
and $\varepsilon_{2}=0.9$. Both probabilities tend to the fixed points
$p_{1}^{*}=0.8$ and $p_{2}^{*}=0.723$.

4. When the herding effect is absent, $p_{1}(t)$ monotonically tends to a fixed point, but $p_{2}(t)$ exhibits an unusual behaviour. At the beginning, for sufficiently long time, $p_{2}(t)$ behaves smoothly. Then suddenly, it starts widely oscillating and continues oscillating for ever. Slightly increasing the herding parameters shifts the beginning of oscillations to larger times, makes the oscillation amplitude smaller, and makes the function $p_{1}(t)$ to also experience everlasting oscillations. The following increase of the herding parameters eliminates oscillations and leads to smooth tendency of both probabilities to different fixed points. The further growth of the herding parameters forces both functions to converge to a consensual fixed point. See Fig. 4.

Figure 4: Probabilities $p_{1}(t)$ (solid line) and $p_{2}(t)$ (dash-dotted
line) for the initial conditions $f_{1}=0.3$, $f_{2}=0$, $q_{1}=0.699$, and
$q_{2}=0.98$. (a) $\varepsilon_{1}=\varepsilon_{2}=0$. The function $p_{1}(t)$
tends to the fixed point $p_{1}^{*}=0.3$, while $p_{2}(t)$, after sufficiently
long time of smooth behaviour, exhibits everlasting oscillations; (b)
$\varepsilon_{1}=0.05$ and $\varepsilon_{2}=0$. Both functions $p_{1}(t)$ and
$p_{2}(t)$ always oscillate; (c) $\varepsilon_{1}=0.05$ and
$\varepsilon_{2}=0.1$. Both functions tend to different fixed points
$p_{1}^{*}=0.322$ and $p_{2}^{*}=0.692$; (d)
$\varepsilon_{1}=\varepsilon_{2}=0.1$. Both probabilities tend to the
consensual fixed point $p_{1}^{*}=p_{2}^{*}=0.98$; (e) $\varepsilon_{1}=0.5$
and $\varepsilon_{2}=0.3$; Both probabilities $p_{1}(t)$ and $p_{2}(t)$
monotonically tend to the consensual fixed point $p_{1}^{*}=p_{2}^{*}=0.98$.
This is an example showing the spontaneous appearance of self-excited
oscillations that can be suppressed by the herding effect.

Figures 5 to 10 illustrate the case of large herding parameters, such that
$\varepsilon_{1}+\varepsilon_{2}>1$. We start with the values
$\varepsilon_{1}=\varepsilon_{2}=1$, so that
$\varepsilon_{1}+\varepsilon_{2}=2$, when the dynamics is described by Eqs.
(36) and then consider diminishing parameters approaching the line
$\varepsilon_{1}+\varepsilon_{2}=1$, where Eq. (35) is valid. The line
$\varepsilon_{1}+\varepsilon_{2}=1$ corresponds to the case of the strongest
herding behavior.

5. For $\varepsilon_{1}+\varepsilon_{2}=2$, the probabilities tend to different fixed points. However, diminishing the herding parameters to the line $\varepsilon_{1}+\varepsilon_{2}=1$ results in the convergence of both functions $p_{j}(t)$ to the common fixed point, as is seen from Fig. 5.

Figure 5: Probabilities $p_{1}(t)$ (solid line) and $p_{2}(t)$ (dash-dotted
line) for the initial conditions $f_{1}=1$, $f_{2}=0.2$, $q_{1}=-0.9$, and
$q_{2}=0.6$. (a) $\varepsilon_{1}=\varepsilon_{2}=1$. The fixed points are
$p_{1}^{*}=0.5$ and $p_{2}^{*}=1$; (b) $\varepsilon_{1}=0.9$ and
$\varepsilon_{2}=0.7$. The consensual fixed point is
$p_{1}^{*}=p_{2}^{*}=0.8$; (c) $\varepsilon_{1}=0.6$ and
$\varepsilon_{2}=0.9$. The consensual fixed point is
$p_{1}^{*}=p_{2}^{*}=0.8$; (d) $\varepsilon_{1}=0.5$ and
$\varepsilon_{2}=0.7$. Both functions tend to $p_{1}^{*}=p_{2}^{*}=0.8$.

6. When $\varepsilon_{1}=\varepsilon_{2}=1$, the probability $p_{1}(t)$ for the agents with long-term memory oscillates from the beginning, while $p_{2}(t)$ tends monotonically to a fixed point. This is contrary to Fig. 3, which is in agreement with Eqs. (36). Diminishing the herding parameters, first, makes both functions permanently oscillating, but then suppresses the oscillations, so that close to the line $\varepsilon_{1}+\varepsilon_{2}=1$ both probabilities tend to fixed points, as is shown in Fig. 6.

Figure 6: Probabilities $p_{1}(t)$ (solid line) and $p_{2}(t)$ (dash-dotted
line) for the initial conditions $f_{1}=0.6$, $f_{2}=1$, $q_{1}=0.39$, and
$q_{2}=-0.9$. (a) $\varepsilon_{1}=\varepsilon_{2}=1$. The function $p_{1}(t)$
oscillates with a constant amplitude and $p_{2}(t)$ tends monotonically to the
fixed point $p_{2}^{*}=0.6$; (b) $\varepsilon_{1}=1$ and
$\varepsilon_{2}=0.9$. Both functions $p_{1}(t)$ and $p_{2}(t)$ permanently
oscillate with constant amplitudes; (c) $\varepsilon_{1}=0.9$ and
$\varepsilon_{2}=0.8$. The functions tend to the fixed points
$p_{1}^{*}=0.265$ and $p_{2}^{*}=0.525$, respectively; (d)
$\varepsilon_{1}=0.3$ and $\varepsilon_{2}=0.9$. The fixed points are
$p_{1}^{*}=0.455$ and $p_{2}^{*}=0.552$.

7. For the case $\varepsilon_{1}=\varepsilon_{2}=1$, the probability $p_{1}(t)$ for the agents with long-term memory permanently oscillates, at the beginning a little, but then the oscillation amplitude increases, while $p_{2}(t)$ monotonically tends to a fixed point. Diminishing the herding parameters leads to the oscillation of both probabilities and shifts the start of oscillations to larger times. Then, reducing the herding parameters to the line $\varepsilon_{1}+\varepsilon_{2}=1$ suppresses the oscillations and forces both functions to converge, first to different fixed points and then to the same consensual fixed point, which is illustrated in Fig. 7.

Figure 7: Probabilities $p_{1}(t)$ (solid line) and $p_{2}(t)$ (dash-dotted
line) for the initial conditions $f_{1}=0.1$, $f_{2}=0$, $q_{1}=0.899$, and
$q_{2}=0.93$. (a) $\varepsilon_{1}=\varepsilon_{2}=1$. The function $p_{1}(t)$
oscillates, while $p_{2}(t)$ tends monotonically to the fixed point
$p_{2}^{*}=0.1$; (b) $\varepsilon_{1}=1$ and $\varepsilon_{2}=0.98$. Both
functions $p_{1}(t)$ and $p_{2}(t)$ oscillate; (c) $\varepsilon_{1}=1$ and
$\varepsilon_{2}=0.9$. The probabilities tend to the fixed points
$p_{1}^{*}=0.6$ and $p_{2}^{*}=0.15$; (d) $\varepsilon_{1}=1$ and
$\varepsilon_{2}=0.8$. Both functions converge to the same fixed point
$p_{1}^{*}=p_{2}^{*}=0.93$; (e) $\varepsilon_{1}=0.5$ and
$\varepsilon_{1}=0.8$. The functions converge monotonically to the consensual
fixed point $p_{1}^{*}=p_{2}^{*}=0.93$ .

8. For the herding parameters on or close to the point $\varepsilon_{1}+\varepsilon_{2}=2$, there appear large chaotic fluctuations after rather long time from the beginning of the process. First, only the probability $p_{1}(t)$ starts chaotically oscillating, while $p_{2}(t)$ smoothly tends to a fixed point. Diminishing the herding parameters results in the sudden chaotic oscillation of both functions. But the further decrease of the herding parameters approaching the line $\varepsilon_{1}+\varepsilon_{2}=1$ , first, makes the oscillations periodic and then eliminates oscillations at al., so that both functions converge to the consensual fixed point (see Fig. 8).

Figure 8: Probabilities $p_{1}(t)$ (solid line) and $p_{2}(t)$ (dash-dotted
line) for the initial conditions $f_{1}=0.3$, $f_{2}=0$, $q_{1}=0.699$, and
$q_{2}=0.99$. (a) $\varepsilon_{1}=\varepsilon_{2}=1$. The function
$p_{1}(t)$, after quite long time of smooth behaviour, suddenly starts
chaotically oscillating, but $p_{2}(t)$ goes monotonically to the fixed point
$p_{2}^{*}=0.3$; (b) $\varepsilon_{1}=1$ and $\varepsilon_{2}=0.95$. Both
functions $p_{1}(t)$ and $p_{2}(t)$, after rather long time of smooth
behaviour, start chaotically oscillating; (c) $\varepsilon_{1}=0.94$ and
$\varepsilon_{2}=0.98$. Both functions $p_{1}(t)$ and $p_{2}(t)$, after the
initial smooth behaviour, start periodic oscillations; (d)
$\varepsilon_{1}=0.9$ and $\varepsilon_{2}=0.8$. Oscillations are suppressed
and the functions $p_{1}(t)$ and $p_{2}(t)$ converge to the consensual fixed
point $p_{1}^{*}=p_{2}^{*}=0.99$; (e) $\varepsilon_{1}=0.5$ and
$\varepsilon_{2}=0.9$. Both functions monotonically converge to the same fixed
point $p_{1}^{*}=p_{2}^{*}=0.99$.

9. At the point $\varepsilon_{1}+\varepsilon_{2}=2$, the probability for the agents with long-term memory begins chaotically oscillating from the beginning, while that for agent with short-term memory monotonically tends to a fixed point. For smaller herding parameters, first, both functions start chaotically oscillating, then chaotic oscillations transfer to periodic fluctuations, and then both functions become monotonic, tending to fixed points. This behaviour is demonstrated in Fig. 9.

Figure 9: Probabilities $p_{1}(t)$ (solid line) and $p_{2}(t)$ (dash-dotted
line) for the initial conditions $f_{1}=0.6$, $f_{2}=1$, $q_{1}=0.3$, and
$q_{2}=-0.999$. (a) $\varepsilon_{1}=\varepsilon_{2}=1$. The function
$p_{1}(t)$ starts chaotically oscillating from the very beginning, while
$p_{2}(t)$ monotonically tends to the fixed point $p_{2}^{*}=0.6$; (b)
$\varepsilon_{1}=1$ and $\varepsilon_{2}=0.8$. Both functions $p_{1}(t)$ and
$p_{2}(t)$ oscillate chaotically; (c) $\varepsilon_{1}=0.9$ and
$\varepsilon_{2}=0.8$. Both functions $p_{1}(t)$ and $p_{2}(t)$ exhibit
periodic oscillations; (d) $\varepsilon_{1}=0.9$ and $\varepsilon_{2}=0.7$.
After the initial period of oscillations, both functions tend to the fixed
points $p_{1}^{*}=0.204$ and $p_{2}^{*}=0.468$; (e) $\varepsilon_{1}=0.2$ and
$\varepsilon_{2}=1$. Both functions monotonically converge to the fixed points
$p_{1}^{*}=0.485$ and $p_{2}^{*}=0.6$.

10. One more example of chaotic fluctuations for large herding parameters, which, first, become periodic and then are suppressed when moving from the point $\varepsilon_{1}+\varepsilon_{2}=2$ to the line $\varepsilon_{1}+\varepsilon_{2}=1$. The transformation of dynamic regimes, when reducing the herding parameters, is shown in Fig. 10.

Figure 10: Probabilities $p_{1}(t)$ (solid line) and $p_{2}(t)$ (dash-dotted
line) for the initial conditions $f_{1}=0.2$, $f_{2}=0$, $q_{1}=-0.1$, and
$q_{2}=0.999$. (a) $\varepsilon_{1}=\varepsilon_{2}=1$. The function
$p_{1}(t)$ experiences chaotic oscillations, while $p_{2}(t)$ monotonically
tends to the fixed point $p_{2}^{*}=0.2$; (b) $\varepsilon_{1}=1$ and
$\varepsilon_{2}=0.7$. Both functions $p_{1}(t)$ and $p_{2}(t)$ chaotically
oscillate; (c) $\varepsilon_{1}=\varepsilon_{2}=0.9$. Both functions
$p_{1}(t)$ and $p_{2}(t)$ demonstrate periodic oscillations; (d)
$\varepsilon_{1}=0.9$ and $\varepsilon_{2}=0.8$. After initial fluctuations,
the functions tend to the fixed points $p_{1}^{*}=0.685$ and
$p_{2}^{*}=0.308$; e) $\varepsilon_{1}=0.1$ and $\varepsilon_{2}=1$. Both
functions monotonically tend to the fixed points $p_{1}^{*}=0.278$ and
$p_{2}^{*}=0.2$.

##  5 Conclusion

We developed a model of a society composed of intelligent agents making
decisions by choosing between several alternatives. The choice is based of
three attributes, utility, attraction, and herding. By varying the herding
parameters there appear different regimes of decision making dynamics. A
detailed investigation of possible regimes is done for the time dependence of
probabilities of choosing one of two alternatives. Two groups of decision
makers are considered, one group consisting of agents with long-term memory
and the other group composed of agents with short-term memory.

The most interesting observation is the sudden appearance of oscillations that
can be either periodic or chaotic. They can arise either from the beginning of
the decision process or after rather long time of smooth dynamics. These
suddenly arising waves are self organized, since the society is a closed
system having no external perturbations. The sudden growth of the self-excited
waves is an interesting example of the appearance of periodic and chaotic
dynamics in a closed system. The herding effect smooths the dynamics and even
can suppress the self-excited waves.

The aim of the paper is the analysis of admissible dynamic regimes in a
complex society. Qualitatively similar dynamics regimes can occur in real
social systems. As illustrations, it is possible to adduce several dynamic
trends from Google analogous to the evolution curves in our figures. Thus,
Fig. 11 demonstrates the curve rising to a limiting value, with small random
fluctuations. Figure 12 depicts the decreasing tendency to a limit. In Fig.
13, we see almost periodic motion. And Fig. 14 displays a kind of chaotic
behavior. It looks that all types of behavior exhibited by the model we
considered can happen in realistic social systems.

Figure 11: Interest to USD in Canada over time from 2004 up to now. (Increases
to its limit).

Figure 12: Interest to car prices worldwide over time from 2004 up to now.
(Decreases to its limit).

Figure 13: Interest to football game in USA over time from 2004 up to now.
(Oscillations).

Figure 14: Interest to liberalism as a political philosophy in Science
category over time for 5 year up to now. (Chaotic behavior).

CRediT authorship contribution statement

V.I. Yukalov and E.P. Yukalova equally contributed to the paper: Concept,
Design, Analysis, Writing, or revision of the manuscript.

Declaration of competing interest

The authors declare that they have no known competing financial interests or
personal relationships that could have appeared to influence the work reported
in this paper.

## References

  * [1] T. Parsons, The Social System, Taylor and Francis, London, 2005.
  * [2] M. Perc, J. Gómez-Gardeñes, A. Szolnoki, L.M. Floria, Y. Moreno, Evolutionary dynamics of group interactions on structured populations: a review, J. Roy. Soc. Interface 10 (2013) 20120997.
  * [3] M. Perc, J.J. Jordan, D.G. Rand, Z. Wang, S. Boccaletti, A. Szolnoki, Statistical physics of human cooperation, Phys. Rep. 687 (2017) 1–51.
  * [4] M. Jusup, P. Holme, K. Kanazawa, M. Takayasu, I. Romic, Z. Wang, S. Gecek, T. Lipic, B. Podobnik, L. Wang, W. Luo, T. Klanjscek, J. Fan, S. Boccaletti, M. Perc, Social physics, arXiv:2110.01866 (2021).
  * [5] J.P. Scott, Social Network Analysis, Sage, Thousand Oaks, 2000.
  * [6] R. Albert, A.L. Barabasi, Statistical mechanics of complex networks, Rev. Mod. Phys. 74 (2002) 47–98.
  * [7] S. Boccaletti, V. Latora, Y. Moreno, M. Chavez, D.U. Hwanga, Complex networks: Structure and dynamics, Phys. Rep. 424 (2006) 175–308.
  * [8] R. Van Meter, Quantum Networking, Hoboken, Wiley, 2014.
  * [9] D.Y. Kenett, M. Perc, S. Boccaletti, Networks of networks: An introduction, Chaos Solit. Fract. 80 (2015) 1–6.
  * [10] A.S. da Mata, Complex networks: A mini review, Braz. J. Phys. 50 (2020) 658–672.
  * [11] N. Nilsson, Artificial Intelligence: A New Synthesis, Morgan Kaufmann, San Francisco, 1998.
  * [12] D. Poole, A. Mackworth, R. Goebel, Computational Intelligence: A Logical Approach, Oxford University, New York, 1998.
  * [13] G.F. Luger, W.A. Stubblefield, Artificial Intelligence: Structures and Strategies for Complex Problem Solving, Benjamin Cummings, Redwood City, 2004.
  * [14] E. Rich, K. Knight, S.B. Nair, Artificial Intelligence, McGraw Hill, New Delhi, 2009.
  * [15] S.J. Russell, P. Norvig, Artificial Intelligence: A Modern Approach, Hoboken, Pearson, 2021.
  * [16] V.I. Yukalov, D. Sornette, Quantum decision theory as quantum theory of measurement, Phys. Lett. A 372 (2008) 6867–6871.
  * [17] V.I. Yukalov, D. Sornette, Scheme of thinking quantum systems, Laser Phys. Lett. 6 (2009) 833–839.
  * [18] V.I. Yukalov, D. Sornette, Physics of risk and uncertainty in quantum decision making, Eur. Phys. J. B 71 (2009) 533–548.
  * [19] V.I. Yukalov, D. Sornette, Decision theory with prospect interference and entanglement, Theory Decis. 70 (2011) 283–328.
  * [20] V.I. Yukalov, D. Sornette, Quantitative predictions in quantum decision theory, IEEE Trans. Syst. Man Cybern. Syst. 48 (2018) 366–381.
  * [21] V.I. Yukalov, Evolutionary processes in quantum decision theory, Entropy 22 (2020) 681.
  * [22] V.I. Yukalov, Tossing quantum coins and dice, Laser Phys. 31 (2021) 055201.
  * [23] V.I. Yukalov, E.P. Yukalova, D. Sornette, Role of collective information in networks of quantum operating agents, Physica A 598 (2022) 127365.
  * [24] V.I. Yukalov, Quantification of emotions in decision making, Soft Comput. 26 (2021) 2419–2436.
  * [25] V.I. Yukalov, A resolution of St. Petersburg paradox, J. Math. Econ. 97 (2021) 102537.
  * [26] P.C. Fishburn, Decision and Value Theory, Wiley, New York, 1964.
  * [27] R.L. Keeney, H. Raiffa, Decisions with Multiple Objectives: Preferences and Value Tradeoffs, Wiley, New York, 1976.
  * [28] S. Frederick, G. Loewenstein, T. O’Donoghue, Time discounting and time preference: A critical review, J. Econ. Liter. 40 (2002) 351–401.
  * [29] V.I. Yukalov, E.P. Yukalova, D. Sornette, Information processing by networks of quantum decision makers, Physica A 492 (2018) 747–766.
  * [30] V.I. Yukalov, Equilibration of quasi-isolated quantum systems, Phys. Lett. A 376 (2012) 550–554.
  * [31] V.I. Yukalov, Decoherence and equilibration under nondestructive measurements, Ann. Phys. (N.Y.) 327 (2012) 253–263.
  * [32] S. Kullback, R.A. Leibler, On information and sufficiency, Ann. Math. Stat. 22 (1951) 79–86.
  * [33] S. Kullback, Information Theory and Statistics, Wiley, New York, 1959.
  * [34] E.D. Martin, The Behavior of Crowds: A Psychological Study, Harper, New York, 1920.
  * [35] M. Sherif, The Psychology of Social Norms, Harper Collins, New York, 1936.
  * [36] N.J. Smelser, Theory of Collective Behavior, Free Press, Glencoe, 1963.
  * [37] R.K. Merton, Social Theory and Social Structure, Free Press, New York, 1968.
  * [38] R.H. Turner, L.M. Killian, Collective Behavior, Prentice-Hall Englewood Cliffs, 1993.
  * [39] E. Hatfield, J.T. Cacioppo, R.L. Rapson, Emotional Contagion, Cambridge University, New York, 1994.
  * [40] M.K. Brunnermeier, Asset Pricing under Asymmetric Information: Bubbles, Crashes, Technical Analysis, and Herding, Oxford University, Oxford, 2001.
  * [41] W.H. Sandholm, Population Games and Evolutionary Dynamics, Massachusetts Institute of Technology, Cambridge, 2010.


# Core packing of well-defined x-ray and NMR structures is the same

Alex T. Grigas  Graduate Program in Computational Biology and Bioinformatics,
Yale University, New Haven, Connecticut, 06520, USA  Integrated Graduate
Program in Physical and Engineering Biology, Yale University, New Haven,
Connecticut, 06520, USA  Zhuoyi Liu  Integrated Graduate Program in Physical
and Engineering Biology, Yale University, New Haven, Connecticut, 06520, USA
Department of Mechanical Engineering and Materials Science, Yale University,
New Haven, Connecticut 06520, USA  Lynne Regan  Institute of Quantitative
Biology, Biochemistry and Biotechnology, Centre for Synthetic and Systems
Biology, School of Biological Sciences, University of Edinburgh  Corey S.
O’Hern  Graduate Program in Computational Biology and Bioinformatics, Yale
University, New Haven, Connecticut, 06520, USA  Integrated Graduate Program in
Physical and Engineering Biology, Yale University, New Haven, Connecticut,
06520, USA  Department of Mechanical Engineering and Materials Science, Yale
University, New Haven, Connecticut 06520, USA  Department of Physics, Yale
University, New Haven, Connecticut 06520, USA  Department of Applied Physics,
Yale University, New Haven, Connecticut 06520, USA

Manuscript Pages: 12

Total Manuscript Figures: 6

Total Manuscript Tables: 0

Supporting Information Pages: 12

Total Supporting Information Figures/Tables: 7/1

Abstract: Numerous studies have investigated the differences and similarities
between protein structures determined by solution NMR spectroscopy and those
determined by x-ray crystallography. A fundamental question is whether any
observed differences are due to differing methodologies, or to differences in
the behavior of proteins in solution versus in the crystalline state. Here, we
compare the properties of the hydrophobic cores of high-resolution protein
crystal structures and those in NMR structures, determined using increasing
numbers and types of restraints. Prior studies have reported that many NMR
structures have denser cores compared to those of high-resolution x-ray
crystal structures. Our current work investigates this result in more detail,
and finds that these NMR structures tend to violate basic features of protein
stereochemistry, such as small non-bonded atomic overlaps and few Ramachandran
and side chain dihedral angle outliers. We find that NMR structures solved
with more restraints, and which do not significantly violate stereochemistry,
have hydrophobic cores that have a similar size and packing fraction as their
counterparts determined by x-ray crystallography at high-resolution. These
results lead us to conclude that, at least regarding the core packing
properties, high-quality structures determined by NMR and x-ray
crystallography are the same, and the differences reported earlier are most
likely a consequence of methodology, rather than fundamental differences
between the protein in the two different environments.

Significance: Dense packing of hydrophobic residues is key to protein
structure and stability. Previously, it has been noted that structures solved
by NMR spectroscopy have denser cores than x-ray crystal structures. Here, we
calculate the core size and packing fraction of NMR structures with
experimental restraints in the PDB. Their cores are typically smaller, but
denser. However, NMR structures with accurate stereochemistry possess core
packing properties that are nearly identical to high-resolution x-ray crystal
structures.

Availability: https://github.com/agrigas115/NMR_code

Keywords: hydrophobic core $|$ protein structure $|$ protein design

##  1 Introduction

We seek to further clarify the key physical properties of proteins that
determine their structure and function. This physics-based approach is
complementary to that of the remarkably successful recent machine learning
methods that predict protein structure from sequence 1, 2. We focus our
analyses on protein hydrophobic cores, because the core region essentially
defines the folded structure and thermodynamic stability of proteins 3, 4.

Over the past three decades, protein structures have been determined primarily
by x-ray crystallography and solution NMR spectroscopy. Prior studies have
reported differences between structures determined by NMR and those determined
by x-ray crystallography 5, 6, 7, 8, 9, 10, 11, 12. A fundamental question is
whether the observed differences are due to differing methodologies, or to
differences in the behavior of proteins in solution versus in the crystalline
state.

For structures determined by x-ray crystallography, it is universally accepted
that a good structure has a low R-factor and that high-resolution structures
represent a more accurate picture of the structure than low-resolution
structures. In previous work, we analyzed a non-redundant set of high-
resolution x-ray crystal structures (with resolution $<1.8$ Å) and found that
the cores (residues with zero solvent accessible surface area) of these
proteins represent about $8\%$ of the total number of residues in the protein.
In addition, the packing fraction, $\phi$, in protein cores is $\sim 0.55\pm
0.01$ 13, 14. In earlier work, using a subset of NMR structures with a large
number of distance restraints per residue, we concluded that the cores of NMR-
determined protein structures had a higher packing fraction than those of
high-resolution protein structures determined by x-ray crystallography 12.
Denser core packing in NMR structures has also been previously reported in
other studies 15, 5. However, for structures determined by NMR, there is no
equivalent universal definition of resolution or metric of quality 16. Here,
we filter the NMR dataset by a wide array of structural metrics and track how
the physical properties of the protein core vary.

Currently, there are approximately $12$,$000$ protein structures solved via
solution-NMR spectroscopy available in the Protein Data Bank (PDB). In this
article, we focus on those structures that have restraints accessible via the
PDB, which is roughly half. We first analyzed the core packing of all
individual structures in these NMR bundles. We found that on average there is
a lower fraction of residues in the core, $f_{c}$ , than in high-resolution
x-ray crystal structures. However, the packing fraction, $\phi$, of the cores
of protein structures determined by NMR is higher than that of the cores
determined by x-ray crystallography. Moreover, these results also hold when
comparing the structures of the same protein, determined by both NMR and x-ray
crystallography.

To investigate the possible origins of the differences seen in NMR core
packing compared to high-resolution x-ray crystal structures, we filtered the
dataset by a variety of structural metrics. We find that NMR bundles solved
with over $10$ distance restraints per residue have a fraction of residues in
the core that is similar to that of high-resolution x-ray crystal structures,
but the average packing fraction of those cores is higher than that observed
in the cores of protein structures determined by x-ray crystallography. We
investigated this phenomenon further by selecting NMR structures that satisfy
common stereochemical metrics, such as low Clashscores and few Ramachandran
and sidechain dihedral angle outliers. Applying these filters results in only
$72$ bundles of the original $6$,$449$. The $f_{c}$ and packing fraction of
the remaining $72$ bundles is the same as that observed in high-resolution
protein x-ray crystal structures. These results suggest that protein
structures are fundamentally the same in solution as in crystals, and that
previously noted differences in the core are a consequence of the structure-
determination methodology.

##  2 Results

First, we assembled datasets of x-ray crystal and NMR protein structures from
the PDB. For the x-ray crystal structures, we used a dataset of $5$,$261$
high-resolution x-ray crystal structures culled from the PDB using PISCES 17,
18 with resolution $<1.8$ Å, sequence identity cutoff of $<20\%$ and R-factor
cutoff of $<0.25$. For the NMR structures, we retrieved all NMR bundles from
the PDB that have available restraint data (i.e. a subset of distance
restraints derived from the nuclear Overhauser effect (NOE), dihedral angle
restraints from J-couplings, and bond-vector orientational restraints from
residual dipolar couplings (RDCs)) using nmr2gmx 19, resulting in $6$,$449$
NMR bundles with a total of $127$,$959$ individual structures. We refer to
these $127$,$959$ NMR structures as the entire NMR dataset. We subsequently
filter this dataset according to the number of NOE distance restraints per
residue and by several stereochemical validation metrics. Our goal is to
understand how the core properties depend on these parameters.

We define core residues as those whose relative solvent accessible surface
area, rSASA, is $\leq 10^{-3}$. We calculate the packing fraction, $\phi$ , of
core residues, which is defined as the ratio of the volume of the residue to
the volume of the Voronoi polyhedron that encloses it 20. We have shown in
previous studies of high-resolution x-ray crystal structures that $\phi$ for
core residues is $\sim 0.55\pm 0.01$. (We also found that $\phi$ increases
with decreasing rSASA 21, until it reaches a plateau at $\phi\sim 0.55$ for
rSASA $<10^{-3}$.) (See Methods for more details.)

In Fig. 1 (A), we compare the distributions of the average packing fraction
$\langle\phi\rangle$ for core residues in the high-resolution x-ray crystal
structure dataset and the entire NMR dataset. For the x-ray crystal
structures, we average over the core residues in each protein structure,
resulting in a single value of $\langle\phi\rangle$ per protein. For the NMR
structures, $\langle\phi\rangle$ is calculated by averaging over all core
residues across all structures in each bundle, resulting in a single
$\langle\phi\rangle$ per protein. We find that for the entire NMR dataset the
average value of $\langle\phi\rangle$ is $0.57$, (two standard deviations
above the value of $0.55$ for x-ray crystal structures), and it is clear in
Fig. 1 (A) that the distribution of $\langle\phi\rangle$ is very broad, and
distinctly skewed towards more densely packed structures.

The average packing fraction in the core, $\langle\phi\rangle$, is essentially
independent of the total number of residues in the protein, $N$. In contrast,
the fraction of the protein that is core, $f_{c}$, increases with increasing
$N$. (See Supporting Information.) Because $f_{c}$ depends on $N$, it is
important to compare datasets with the same length distribution. The proteins
in the NMR data set tend to be smaller than those in the dataset of x-ray
crystal structures, therefore to calculate $P(f_{c})$, we resampled the x-ray
crystal dataset so that it possesses a $P(N)$ that is similar to that for the
NMR dataset for $N<300$. The resampling was obtained by averaging over
$10^{3}$ independent trials. $P(f_{c})$ for the NMR dataset and the resampled
x-ray crystal structure dataset with $N<300$ is plotted in Fig. 1 (B). It is
clear that $P(f_{c})$ has a strong peak for $f_{c}\sim 0$ (i.e. no core
residues) for the NMR dataset. In contrast, $P(f_{c})$ shows that most x-ray
crystal structures of these lengths have a core of $f_{c}\sim 0.05$ on
average. In summary, when we analyze all structures in the PDB determined by
NMR, for which restraint data is available, they on average possess smaller,
more densely packed cores compared to high-resolution x-ray crystal
structures. Interestingly, smaller, but more densely packed cores have also
been found in many computationally generated low-scoring structures with
incorrect backbone placement submitted to the Critical Assessment of protein
Structure Prediction (CASP) competitions 14.

Figure 1: (A) The distribution of the packing fraction of core residues
$P(\langle\phi\rangle)$ for high-resolution x-ray crystal structures (grey
squares) and NMR structures (black circles). For the x-ray crystal structures,
$\langle\phi\rangle$ is averaged over all core residues for each protein and
for the NMR structures, $\langle\phi\rangle$ is averaged over all core
residues for all structures in the bundle for each protein. (B) The
distribution $P(f_{c})$ for the resampled x-ray crystal structures (grey
squares with shading based on the standard deviation) averaged over $10^{3}$
resampling trials and for NMR structures (black circles), where the protein
size is limited to $50<N<300$.

To investigate the discrepancies in the core packing properties in greater
detail, we compared protein structures that had been determined by both x-ray
crystallography and NMR. We assembled a dataset of $702$ pairs with both a
high-resolution x-ray crystal structure ($<2.0$ Å) and an NMR structure with
available restraints, where the pairs have $>90\%$ sequence similarity. (See
Methods.) In Fig. 2, we show a scatter plot of the differences in the average
core packing fraction and fraction of core residues between the NMR and x-ray
crystal structures, where $\Delta\langle\phi\rangle=\langle\phi\rangle_{\rm
NMR}-\langle\phi\rangle_{\rm x-ray}$ and $\Delta\langle f_{c}\rangle=\langle
f_{c}\rangle_{\rm{NMR}}-\langle f_{c}\rangle_{\rm x-ray}$. Consistent with
Fig. 1, the majority of NMR structures have smaller, more densely packed cores
compared to their x-ray crystal structure counterparts, i.e. most data points
occur in the upper left corner. However, there is also data in the three other
quadrants.

Figure 2: Scatter plot of the differences in the packing fraction
$\Delta\langle\phi\rangle$ and fraction of core resides $\Delta f_{c}$ between
the NMR and x-ray crystal structures in the paired dataset. Figure 3: (A) The
distribution $P(N_{\rm{r}}/N)$ of the number of distance restraints per
residue $N_{\rm{r}}/N$ for the NMR dataset. (B) The distribution
$P(\langle\phi\rangle)$ of the average core packing fraction and (C)
distribution $P(f_{c})$ of the fraction of core residues for NMR structures
with $N_{\rm{r}}/N$ greater than a given value (indicated by the color from
$0$ in violet to $30$ in dark red). In (B) and (C), we also show
$P(\langle\phi\rangle)$ and $P(f_{c})$ for the resampled high-resolution x-ray
crystal dataset (grey squares with shading that indicates one standard
deviation).

These results raise the interesting question of whether the differences in the
size and packing fraction of the core represent true differences between
protein structures when they are in crystalline states versus in solution? Or
are the differences a consequence of the methods used to determine the atomic
structure? Since NMR bundles are determined by a set of restraints (i.e.
distance restraints derived from NOEs, dihedral angle restraints from
J-couplings, and bond-vector orientational restraints from RDCs), it is
reasonable to assume that the more restraints used, the higher the quality of
the NMR structure. Note that while requiring a large number of restraints is
necessary, previous studies have shown that it is not sufficient for
identifying high-quality NMR structures 22, 23, 24. Nevertheless, it is of
interest to track how the core properties vary with the number of NOE
measurements per residue. In a typical NMR structure there are many more NOE
distance restraints than any other type of restraint, so we chose to analyse
their effect on core properties.

We show the distribution $P(N_{\rm{r}}/N)$ of the number of NOE distance
restraints per residue in the NMR dataset in Fig. 3 (A). There are many NMR
structures with $N_{\rm{r}}/N\sim 1$, but the vast majority have
$N_{\rm{r}}/N\sim 10$-$25$, and a few have $N_{\rm{r}}/N\gtrsim 30$. In Fig. 3
(B), we plot the distribution of the average packing fraction of core
residues, $P(\langle\phi\rangle)$, in NMR structures with different values of
$N_{\rm{r}}/N$. We find that $P(\langle\phi\rangle)$ does not depend strongly
on $N_{\rm{r}}/N$ and differs from $P(\langle\phi\rangle)$ for high-resolution
x-ray crystal structures for all values of $N_{\rm{r}}/N$. In contrast, the
distribution of the average fraction of core residues, $P(f_{c})$, in NMR
structures becomes increasingly similar to that of the high-resolution x-ray
crystal structure dataset for $N_{\rm{r}}/N\gtrsim 20$, as shown in Fig. 3
(C). There is no correlation between $N$ and $N_{\rm{r}}/N$ and therefore we
used the same resampling of the x-ray crystal structure database in Fig. 1.
The number of NOE distance restraints per residue has been reported to be only
an approximate proxy of the quality of the protein structure, since the NOE
distance restraints can be redundant 25. However, when we perform the same
analyses for the more informative inter-residue distance restraints involving
side chain atoms, the results are the same. (See the Supporting Information.)

We were interested in investigating if the packing fraction of core residues
in NMR structures changed over time as the methodologies for NMR structure
determination have improved 26, 27, 28. In Fig. 4 (A), we show the average
packing fraction of core residues $\langle\phi\rangle$ versus the year that
the NMR structures were deposited in the PDB. We find that
$\langle\phi\rangle$ was largest before $2004$ ($\langle\phi\rangle\sim
0.59$), decreased to $\langle\phi\rangle\sim 0.57$ after 2004, and has
remained nearly constant over the past decade. One possible cause of the
elevated values of $\langle\phi\rangle$ in the NMR structures, relative to the
value of $0.55$ observed in high-resolution x-ray crystal structures, is
overlap of non-bonded atoms, which can increase the apparent packing fraction.
It is therefore of interest to note that the Clashscore, defined by the
MolProbity software 28 as the number of non-bonded atomic overlaps that are
greater than $0.4$ Å per $1$,$000$ atoms, has decreased over time for NMR
structures, with a rapid decrease around $2004$ (the year MolProbity was first
published) and then a slow, but steady decrease towards the average value for
high-resolution x-ray crystal structures over the last decade. (See Fig. 4
(B).) The wwPDB NMR Validation Task Force has emphasized that validation
metrics, such as the Clashscore and Ramachandran and side chain dihedral angle
outliers, used to determine the quality of x-ray crystal structures are
appropriate for assessing NMR structures 16. In particular, these metrics
should be similar for all high-quality protein structures, i.e. structures
that are obtained from x-ray crystallography and solution NMR spectroscopy 29,
26. When we remove NMR structures that have a Clashscore larger than the
average value plus one standard deviation of the Clashscore for high-
resolution x-ray crystal structures (i.e. Clashscore $>10.5$), we find that
$\langle\phi\rangle$ decreases significantly, but $\langle\phi\rangle$ for
these NMR structures with minimal or no overlap of non-bonded atoms, remains
above the value for high-resolution x-ray crystal structures.

Figure 4: (A) The average packing fraction $\langle\phi\rangle$ of core
residues and (B) Clashscore in NMR structures versus the year that they were
deposited in the PDB for the entire NMR dataset (black circles) and those with
a Clashscore that is lower than the average value for high-resolution x-ray
crystal structures plus one standard deviation (grey squares). The average
values and standard deviations of $\langle\phi\rangle$ and Clashscore for
high-resolution x-ray crystal structures are indicated by red dashed lines and
blue shading, respectively.

These results indicate that while removing NMR structures that contain
unphysical atomic overlaps accounts for some of the overpacking found in NMR
structures, $\langle\phi\rangle$ for NMR structures is still above the value
for x-ray crystal structures. Thus, there is another reason, in addition to
atomic overlaps, that is responsible for the elevated packing fraction
observed in some NMR structures. Although many of the metrics are inter-
related, we considered several additional quality metrics, and analysed how
the properties of the cores of the protein structures depend on them. To this
end, we explicitly consider the fraction of Ramachandran and sidechain
dihedral angle outliers, bond length and bond angle outliers, fraction of
backbone chemical shifts assigned, and strength of the magnet used. If we
filter the NMR dataset using each of these metrics one at a time,
$\langle\phi\rangle$ for NMR structures decreases toward the value for high-
resolution x-ray crystal structures and the $P(f_{c})$ distribution for NMR
structures more closely resembles that of the x-ray crystal structure
distribution. However, whether or not the structure was solved with RDCs makes
little difference for the core packing properties. (See Supporting
Information.)

Figure 5: (A) The distribution $P(\langle\phi\rangle)$ of the average packing
fraction for core residues and (B) the distribution $P(f_{c})$ of the fraction
of core residues for high-resolution x-ray crystal structures (grey squares),
the filtered NMR structures (red circles), and the entire NMR dataset (black
circles dashed).

For determining the validation metric cutoffs, we sought a compromise between
strict cutoffs near the high-resolution x-ray crystal structure values and
loose cutoffs that admit more data. When we filter the NMR dataset so that the
NMR structures 1) possess values for the Clashscore, Ramachandran and
sidechain dihedral angle outliers that are within one standard deviation of
the average for the high-resolution x-ray crystal structures (i.e. Clashscore
$<10.5$, Ramachandran dihedral angle outliers $<0.5\%$ and side chain dihedral
angle outliers $<3\%$), 2) have at least $10$ NOE distance restraints per
residue, 3) have $60\%$ or more backbone chemical shifts assigned, and 4) were
solved on a $650$ MHz magnet or stronger, $\langle\phi\rangle$ and $P(f_{c})$
for the NMR structures are nearly identical to those for the high-resolution
x-ray crystal structures as shown in Fig. 5.

These filters are quite stringent, as the filtered dataset contains only $72$
NMR structures of the $6$,$449$ from the entire NMR dataset. The average
packing fraction of core residues in this filtered NMR dataset is
$\langle\phi\rangle=0.55\pm 0.02$ compared to the value of
$\langle\phi\rangle=0.55\pm 0.01$ for the high-resolution x-ray crystal
structure dataset. The distribution of the fraction of core residues is also
similar for the filtered NMR structures and high-resolution x-ray crystal
structures. (See Figure 5.) Similar results can be achieved when only the
stereochemical criteria are considered. (See Supporting Information.) None of
the filtered NMR bundles is included the x-ray/NMR paired dataset.

Figure 6: (A) Average packing fraction of core residues $\langle\phi\rangle$
and average fraction of core residues $\langle f_{c}\rangle$ for the x-ray
crystal structure with PDBID: 2CRW (x), the NMR bundle with PDBID: 2CZN (+),
and unrestrained (upward triangle) and restrained (circles) MD simulations
colored by the average fraction of NOE distance restraints $\langle
F_{r}\rangle$ that are satisfied. (B) The average Clashscore plotted versus
the average fraction of NMR distance restraints satisfied $\langle
F_{r}\rangle$ for the systems in (A). The data for the NMR bundle with a high
Clashscore of $124$ is not shown.

To investigate the origin of the deviation in the packing properties between
NMR and x-ray crystal structures, we used GROMACS to carry out all-atom MD
simulations of a globular protein, the chitin binding domain of chitinase from
Pyrococcus furiosus, which has an x-ray crystal structure with a resolution of
$1.7$ Å (PDBID: 2CWR) and an NMR bundle with a large number of distance
restraints (PDBID: 2CZN). There are $38$ individual models in the NMR bundle,
solved with an average of $25$ distance restraints per residue. Nevertheless,
the bundle has a high Clashscore of $124$. We used MD simulations to
investigate if we can remove the atomic overlaps in the NMR bundle and match
the NMR distance restraints to the same extent as the deposited bundle. When
retrieving NOE distance restraints from the PDB, the restraints are reported
as a lower bound, a first upper bound, and a second upper bound on the
pairwise atomic distances. Here, we define a distance restraint as being
satisfied if the distance is less than the first upper bound. In the case of
ambiguous NOE restraints, where the NOE restraint is attributed to more than a
single pair of atoms (e.g. methyl hydrogens), we count the restraint as
satisfied if at least one of the set of pairs has a distance less than the
first upper bound. To assess the resulting structures, we then calculated
$\langle\phi\rangle$ and $\langle f_{c}\rangle$ over the MD simulation
trajectory.

In Fig. 6 (A), we plot the average fraction of core residues $\langle
f_{c}\rangle$ versus the average packing fraction of core residues
$\langle\phi\rangle$ from the MD simulations and for the x-ray crystal and NMR
structures. For this particular protein, the NMR bundle has both a higher
packing fraction and a larger fraction of core residues than the x-ray crystal
structure, but it has a large Clashscore and many Ramachandran and sidechain
dihedral angle outliers. The color scale in Fig. 6 (A) indicates the average
fraction $F_{r}$ of the distance restraints that are satisfied (using the
first upper bound).

We first ran unrestrained MD simulations, starting from a structure in the NMR
bundle. This data point is plotted as a triangle in Fig. 6 (A). In the
unrestrained MD simulations, $85\%$ of the distance restraints are satisfied
compared to $\sim 93\%$ in the x-ray crystal structure (x) and $\sim 100\%$ in
the NMR bundle (+). Note that in the unrestrained MD simulations
$\langle\phi\rangle$ decreases from $\langle\phi\rangle\sim 0.63$ for the
initial NMR structure to the value for the x-ray crystal structure,
$\langle\phi\rangle\sim 0.57$. However, the fraction of core residues,
$\langle f_{c}\rangle\sim 0.03$, is much smaller than that for the initial NMR
bundle with $\langle f_{c}\rangle\sim 0.11$ and the x-ray crystal structure
with $f_{c}\sim 0.08$. The unrestrained MD simulation was run for $500$ ns,
during which time the packing properties rapidly equilibrate, and therefore
the restrained MD simulations described below were run for shorter times. (See
Supporting Information.)

We also ran MD simulations with restraints to recapitulate the NOE distances.
GROMACS enforces NMR distance restraints using the following potential energy:

$V_{dr}(r_{ij})=\begin{cases}\frac{k}{2}(r_{ij}-r_{0})^{2},&r_{ij}<r_{0},\\\
0,&r_{0}\leq r_{ij}<r_{1},\\\ \frac{k}{2}(r_{ij}-r_{1})^{2},&r_{1}\leq
r_{ij}<r_{2},\\\ \frac{k}{2}(r_{2}-r_{1})(2r_{ij}-r_{2}-r_{1}),&r_{2}\leq
r_{ij},\\\ \end{cases}$ (1)

where $r_{ij}$ is the separation between NOE atom pairs $i$ and $j$, $r_{0}$,
$r_{1}$, and $r_{2}$ are the lower bound, first upper bound, and second upper
bound distances for the NOE atom pairs. Additionally, the average separation
$\langle r_{ij}\rangle$ is weighted by a factor of $r^{-3}_{ij}$ multiplied by
an exponentially decaying function of time with decay constant $\tau$. (See
Methods.) This restraint potential is similar to the potentials used in NMR
structure determination software 30, 31, 32. Each restrained simulation was
initialized using a single structure from the NMR bundle and run for $20$ ns
following energy minimization, and the NOE pair separations and packing
properties were obtained every $0.1$ ns and averaged. We considered several
spring constants $k=10^{3}$, $10^{4}$, and $10^{5}$ kJ/mol/nm2 and a decay
time $\tau=0.5$ ns, which increases the average fraction $\langle
F_{r}\rangle$ of distance restraints that are satisfied from $0.86$ to $0.90$.
For $k>10^{5}$ kJ/mol/nm2, $\langle F_{r}\rangle$ saturates and cannot reach
$\langle F_{r}\rangle\sim 1$, which is the value for the NMR bundle deposited
in the PDB. (See Supporting Information.) Restrained MD simulations using the
flat-bottomed potential in Eq. 1 yield cores with similar values of $\langle
f_{c}\rangle$ and $\langle\phi\rangle$ as the x-ray crystal structure. To
obtain $\langle F_{r}\rangle\sim 1$, we implemented a new flat-bottomed
distance restraint potential, where $r_{1}$ is set half way between the
original lower bound and original first upper bound, using spring constants
$k=10^{3}$, $10^{4}$, and $10^{5}$ kJ/mol/nm2 and a decay time $\tau=0.5$ ns
and a final highly restrained simulation with a spring constant $k=10^{5}$
kJ/mol/nm2 and $\tau=0$ ns. (See Methods.) With these modifications to the
restraint potential and $\tau=0$ ns, we are able to obtain structures that
have $\langle F_{r}\rangle\sim 1$ and values of $\langle f_{c}\rangle$ and
$\langle\phi\rangle$ that match those for the NMR bundle. However, in Fig. 6
(B), we show that as $\langle F_{r}\rangle$ increases toward $1$, the
Clashscore rises steeply, reaching values that are indicative of a poor
quality protein structure. In addition, as shown in the Supporting
Information, the same behavior is found for the fraction of Ramachandran and
side chain dihedral angle outliers. The final increase in $\langle
F_{r}\rangle$ toward $1$ is caused by the increase in the number of satisfied
inter-residue distance restraints involving side chain atoms. (See the
Supporting Information.) By enforcing the distance restraints such that
$\langle F_{r}\rangle\sim 1$ in the MD simulations, we obtain the same core
packing properties as the deposited NMR bundle, but these structures violate
metrics of protein stereochemistry. Thus, matching all of the deposited NOE
restraints to the same extent as the deposited NMR bundle is likely the cause
of the overpacking observed in this particular NMR bundle.

##  3 Discussion

An important question in protein science is to what extent are the structures
of proteins in the crystalline state different from those of proteins in
solution. Here, we have shown that most NMR structures deposited in the PDB
possess smaller cores, and that these cores are more densely packed than those
observed in high-resolution x-ray crystal structures. More densely packed
cores in NMR structures have also been reported previously for smaller
datasets 12, 5, 15. However, we find that by only considering NMR bundles that
have a large number of restraints and that pass protein stereochemistry
validation metrics (i.e. those that have small Clashscores and few
Ramachandran and sidechain dihedral angle outliers), both the core size and
packing fraction of these NMR structures match the values found for high-
resolution x-ray crystal structures. The validation metrics are quite
stringent: Only $72$ NMR bundles of the original $6$,$449$ bundles in the
dataset pass the validation metrics. We interpret these results to show that
the packing properties of the cores of proteins in the crystalline state
(determined by high-resolution x-ray crystallography) and the cores of
proteins in solution (determined by state-of-the-art NMR spectroscopy) are the
same. However, since none of the filtered NMR structures have high-resolution
x-ray pairs, it remains possible that particular proteins may have non-zero
$\Delta\langle\phi\rangle$ or $\Delta f_{c}$.

A recent study reported that many NMR structures have surface loops that are
often under-restrained and too floppy 33. This observation is consistent with
our finding that the cores of NMR structures solved with few distance
restraints are typically too small. In particular, inappropriately positioned
surface residues will affect the solvent accessibility of the interior
residues.

Over-packing occurs simultaneously with poor protein stereochemistry (large
non-bonded atomic clashes and many dihedral angle outliers) in NMR structures
in the PDB. However, this result does not rule out the possibility that denser
core packing than that found in high-resolution x-ray crystal structures could
occur in principle without violating protein stereochemistry. For example,
packings of amino acid-shaped particles can become denser with the addition of
thermal fluctuations, in contrast to packings that are rapidly compressed and
energy-minimized 12. Therefore, we carried out MD simulations with NOE
distance restraints to attempt to remove non-bonded atomic overlaps and
dihedral angle outliers and then measured the resulting core packing
properties. Our MD simulations on a particular globular protein, the chitin
binding domain of chitinase from Pyrococcus furiosus, showed that forcing the
NOE distance restraints to be satisfied to the same extent as the deposited
NMR bundle leads to incorrect core packing features, non-bonded atomic
clashes, and dihedral angle outliers.

The NOE restraint data clearly contains important structural information, as
the NMR and the x-ray crystal structures have similar overall folds (with an
average Cα RMSD $=1.6$ Å) and loosely enforcing all of the NOE restraints in
the MD simulations leads to core packing properties that are similar to the
high-resolution x-ray crystal structure. However, enforcing the NOE distance
restraints to the same extent of the deposited NMR bundle appears to be
inconsistent with protein steroechemistry. The fact that simultaneously
satisfying all of the NOE distance restraints leads to violations of protein
stereochemistry (for this particular protein) is likely related to the fact
that the conversion of NOE measurements to distance bounds is approximate and
NOEs report on time-averaged atomic separations, which can be challenging to
implement in protein structure determination 23, 34, 35. We emphasize that our
MD simulations were performed on a single protein whose structure was
determined by both x-ray crystallography and NMR (from quadrant 2 of Fig. 2).
In future studies, it will be important to study protein structure pairs that
represent all four quadrants from Fig. 2 and to determine what restrained MD
protocol would be needed to improve the modeling of the protein core.
Additionally, our unrestrained MD simulations show that the core packing
fraction is recapitulated using the Amber99SB-ILDN forcefield and TIP3P water
model, but the core size is not. Therefore, we will explore whether the same
result is observed with different force fields and solvation models, which
would indicate how well the hydrophobic effect is currently captured in MD
simulations of proteins.

##  4 Materials and Methods

### Datasets

We collected two datasets of protein structures from the Protein Data Bank
(PDB). We first compiled a dataset of $5$,$261$ high-resolution x-ray crystal
structures from the PDB using PISCES 17, 18 with resolution $<1.8$ Å, a
sequence identity cutoff of $<20\%$, an R-factor cutoff of $<0.25$, and
lengths greater than $50$ residues. We also assembled a large dataset of NMR
structures from the PDB. We identified all NMR bundles that have available
restraint data in the NMR-STAR format using nmr2gmx 19, resulting in $6$,$499$
NMR bundles, each with at least $10$ structures in the bundle. In this
dataset, $5$,$830$ bundles have NOE distance restraints, $3$,$994$ have
dihedral angle restraints, and $507$ have RDC restraints. The NMR dataset was
not filtered for similarity in amino acid sequence, and thus it has multiple
bundles for some proteins. PyPDB was used to retrieve all of the validation
metrics for each protein from the PDB 36. In addition, PyPDB was used to
identify a set of NMR and x-ray crystal structure pairs. We selected pairs for
which the x-ray crystal structure sequence matched $>90\%$ via MMSeqs2 37, the
lengths differ by $<10\%$, the x-ray crystal structure resolution is $\leq
2.0$ Å, and the NMR bundle has restraint data. This dataset has $702$ pairs,
made up of $514$ x-ray crystal structures and $525$ NMR structures, as some
NMR structures have multiple x-ray crystal structure matches and vice versa.

### rSASA

To identify core residues, we measured each residue’s solvent accessible
surface area (SASA). To calculate SASA, we use the Naccess software package
38, which implements an algorithm originally proposed by Lee and Richards 39.
To normalize the SASA, we take the ratio of the SASA within the context of the
protein ($\text{SASA}_{\text{context}}$) and the SASA of the same residue (X)
extracted from the protein structure as a dipeptide (Gly-X-Gly) with the same
backbone and side-chain dihedral angles:

$\text{rSASA}=\frac{\text{SASA}_{\text{context}}}{\text{SASA}_{\text{dipeptide}}}.$
(2)

Core residues are classified as those that have ${\rm rSASA}\leq 10^{-3}$. The
choice of this threshold has been discussed in several previous studies 21,
20.

### Packing Fraction

A characteristic measure of the packing in an atomic system is the packing
fraction. The packing fraction of residue $\mu$ is defined by

$\phi_{\mu}=\frac{\nu_{\mu}}{V_{\mu}},$ (3)

where $\nu_{\mu}$ is the non-overlapping volume of residue $\mu$ and $V_{\mu}$
is the volume of the Voronoi cell surrounding residue $\mu$. The volume of the
residue is determined by the atom sizes and their locations. The particular
atom sizes that we used have been previously selected from the range of
literature values to reproduce the backbone and side chain dihedral angle
distributions found in high-resolution x-ray crystal structures 40, 41, 42.
The Voronoi cell represents the local free space around the residue. We
calculate the non-overlapping residue volume with a grid-based volume
estimation. To calculate the Voronoi polyhedra for a protein structure, we use
surface Voronoi tessellation, which defines a Voronoi cell as the region of
space in a given system that is closer to the bounding surface of the residue
than to the bounding surface of any other residue in the system. We calculate
the surface Voronoi tessellations using the Pomelo software package 43. This
software approximates the bounding surfaces of each residue by triangulating
points on the residue surfaces. We find that using $\sim 400$ points per atom,
or $\sim 6400$ surface points per residue, gives an accurate representation of
the Voronoi cells and the results do not change if more surface points are
included.

### Molecular Dynamics Simulations

The unrestrained and restrained all-atom molecular dynamics (MD) simulations
of the globular protein, the chitin binding domain of chitinase from
Pyrococcus furiosus, were performed using the GROMACS software package 44. The
MD simulations used a cubic box with edge lengths $63$ Å, which ensures that
the protein is $>10$ Å away from the box edges, and the box is filled with
water molecules modeled using TIP3P 45 at neutral pH and $0.15$M NaCl.
Periodic boundary conditions are applied in the $x$-, $y$-, and
$z$-directions. Short-range van der Waals and screened Coulomb interactions
were truncated at $1.2$ nm, while longer-ranged electrostatics were tabulated
using the Particle Mesh Ewald summation method. We performed two short energy
minimization runs to relax the protein first and then the water molecules and
protein together using the steepest decent method. We then simulated the
system for $20$ ns in the NPT ensemble at temperature $T=300$K and pressure
$P=1$ bar using the weakly coupled Berendsen thermostat and barostat and
calculated the $\langle\phi\rangle$ and $f_{c}$ every $0.1$ ns. The time
constant of the Berendsen thermostat was set to $2$ ps and the isothermal
compressibility for the Berendsen barostat was set to $4.5\times 10^{-5}\rm
bar^{-1}$. (See the Supporting Information for more details about the
equilibration procedure.) The equations of motion for the atomic coordinates
and velocities were integrated using a leapfrog algorithm with a $1$ fs time
step. For both the unrestrained and restrained MD simulations, we used the
AMBER99SB-ILDN force field 46, 47. For the restrained MD simulations, we
implemented the flat-bottomed spring potential $V_{dr}(r_{ij})$ given in Eq. 1
with spring constants $k=10^{3}$, $10^{4}$, and $10^{5}$ kJ/mol/nm2 to enforce
the NOE distance restraints between atom pairs $i$ and $j$. We calculate the
time average NOE pairwise atomic separations as:

$\langle r_{ij}\rangle=\left\langle\left(r^{-3}_{ij}(t-\Delta t)e^{(-\Delta
t/\tau)}+r^{-3}(t)\left(1-e^{(-\Delta
t/\tau)}\right)\right)^{-1/3}\right\rangle_{t},$ (4)

where $\Delta t$ is the time step and $\tau$ is the decay time. $\tau=0.5$ ns
for all restraint simulations, except for the most highly restrained
simulation. This time-averaging method has been suggested to more accurately
represent NOE distance restraints and allows for larger fluctuations in the
atomic separations if NOE distance restraints are incompatible. To further
increase the fraction of NOE restraints that are satisfied during the MD
simulations to $\langle F_{r}\rangle\sim 1$, we also implemented a narrower
flat-bottomed spring potential, where the first upper bound $r_{1}$ is set to
the separation between the NOE atom pairs of the initial NMR structure, and
$\tau=0$ when calculating $\langle r_{ij}\rangle$.

## References

  * 1 Pereira J, Simpkin AJ, Hartmann MD, Rigden DJ, Keegan RM, Lupas AN (2021) High-accuracy protein structure prediction in CASP14. Proteins Struct Funct Bioinf 89:1687–1699.
  * 2 Jumper J, Evans R, Pritzel A, Green T, Figurnov M, Ronneberger O, Tunyasuvunakool K, Bates R, Žídek A, Potapenko A, Bridgland A, Meyer C, Kohl SAA, Ballard AJ, Cowie A, Romera-Paredes B, Nikolov S, Jain R, Adler J, Back T, Petersen S, Reiman D, Clancy E, Zielinski M, Steinegger M, Pacholska M, Berghammer T, Bodenstein S, Silver D, Vinyals O, Senior AW, Kavukcuoglu K, Kohli P, Hassabis D (2021) Highly accurate protein structure prediction with AlphaFold. Nature 596:583–589.
  * 3 Dill KA (1990) Dominant forces in protein folding. Biochemistry 29:7133–7155.
  * 4 Pace CN, Fu H, Fryar KL, Landua J, Trevino SR, Shirley BA, Hendricks MM, Iimura S, Gajiwala K, Scholtz JM, Grimsley GR (2011) Contribution of hydrophobic interactions to protein stability. J Mol Biol 408(3):514–528.
  * 5 Garbuzynskiy SO, Melnik BS, Lobanov MY, Finkelstein AV, Galzitskaya OV (2005) Comparison of X-ray and NMR structures: Is there a systematic difference in residue contacts between X-ray and NMR-resolved protein structures? Proteins Struct Funct Bioinf 60:139–147.
  * 6 Andrec M, Snyder DA, Zhou Z, Young J, Montelione GT, Levy RM (2007) A large data set comparison of protein structures determined by crystallography and NMR: Statistical test for structural differences and the effect of crystal packing. Proteins Struct Funct Bioinf 69:449–465.
  * 7 Schneider M, Fu X, Keating AE (2009) X-ray vs. NMR structures as templates for computational protein design. Proteins Struct Funct Bioinf 77: 97–110.
  * 8 Sikic K, Tomic S, Carugo O (2010) Systematic comparison of crystal and NMR protein structures deposited in the Protein Data Bank. Open Biochem J 4:83–95.
  * 9 Mao B, Tejero R, Baker D, Montelione GT (2014) Protein NMR structures refined with Rosetta have higher accuracy relative to corresponding X-ray crystal structures. J Am Chem Soc 136:1893–1906.
  * 10 Everett JK, Tejero R, Murthy SBK, Acton TB, Aramini JM, Baran MC, Benach J, Cort JR, Eletsky A, Forouhar F, Guan R, Kuzin AP, Lee HW, Liu G, Mani R, Mao B, Mills JL, Montelione AF, Pederson K, Powers R, Ramelot T, Rossi P, Seetharaman J, Snyder D, Swapna GVT, Vorobiev SM, Wu Y, Xiao R, Yang Y, Arrowsmith CH, Hunt JF, Kennedy MA, Prestegard JH, Szyperski T, Tong L, Montelione GT (2016) A community resource of experimental data for NMR / X-ray crystal structure pairs. Protein Sci, 25:30–45.
  * 11 Leman JK, D’Avino AR, Bhatnagar Y, Gray JJ (2018) Comparison of NMR and crystal structures of membrane proteins and computational refinement to improve model quality. Proteins Struct Funct Bioinf 86:57–74.
  * 12 Mei Z, Treado JD, Grigas AT, Levine ZA, Regan L, O’Hern CS (2020) Analyses of protein cores reveal fundamental differences between solution and crystal structures. Proteins Struct Funct Bioinf 88:1154–1161.
  * 13 Gaines JC, Smith WW, Regan L, O’Hern CS (2016) Random close packing in protein cores. Phys Rev E 93:032415.
  * 14 Grigas AT, Mei Z, Treado JD, Levine ZA, Regan L, O’Hern CS (2020) Using physical features of protein core packing to distinguish real proteins from decoys. Protein Sci 29:1931–1944.
  * 15 Ratnaparkhi GS, Ramachandran S, Udgaonkar JB, Varadarajan R. (1998) Discrepancies between the NMR and X-ray structures of uncomplexed barstar: Analysis suggests that packing densities of protein structures determined by NMR are unreliable. Biochemistry, 37:6958–6966.
  * 16 Montelione GT, Nilges M, Bax Ad, Güntert P, Herrmann T, Richardson JS, Schwieters CD, Vranken WR, Vuister GW, Wishart DS, Berman HM, Kleywegt GJ, Markley JL (2013) Recommendations of the wwPDB NMR Validation Task Force. Structure 21:1563–1570.
  * 17 Wang G, Dunbrack RL (2003) PISCES: A protein sequence culling server. Bioinformatics 19:1589–1591.
  * 18 Wang G, Dunbrack RL (2005) PISCES: Recent improvements to a PDB sequence culling server. Nucleic Acids Res, 33:W94–W98.
  * 19 Sinelnikova A, van der Spoel D (2021) NMR refinement and peptide folding using the GROMACS software. J Biomol NMR, 75:143–149.
  * 20 Treado JD, Mei Z, Regan L, O’Hern CS (2019) Void distributions reveal structural link between jammed packings and protein cores. Phys Rev E, 99:022416.
  * 21 Gaines JC, Acebes C, Virrueta A, Butler M, Regan L, O’Hern CS (2018) Comparing side chain packing in soluble proteins, protein-protein interfaces and transmembrane proteins. Proteins Struct Funct Bioinf 86:581–591.
  * 22 Spronk CAEM, Nabuurs SB, Bonvin AMJJ, Krieger E, Vuister GW, Vriend G (2003) The precision of NMR structure ensembles revisited. J Biomol NMR, 25:225–234.
  * 23 Zagrovic B, van Gunsteren WF. Comparing atomistic simulation data with the NMR experiment: How much can NOEs actually tell us? Proteins Struct Funct Bioinf 63:210–218.
  * 24 Rosato A, Tejero R, Montelione GT (2013) Quality assessment of protein NMR structures. Curr Opin Struct Biol 23:715–724.
  * 25 Nabuurs SB, Spronk CAEM, Krieger E, Maassen H, Vriend G, Vuister GW (2003) Quantitative evaluation of experimental NMR restraints. J Am Chem Soc 125:12026–12034.
  * 26 Vuister GW, Fogh RH, Hendrickx PMS, Doreleijers JF, Gutmanas A (2014) An overview of tools for the validation of protein NMR structures. J Biomol NMR 58:259—285.
  * 27 Buchner L, Güntert P (2015) Increased reliability of nuclear magnetic resonance protein structures by consensus structure bundles. Structure 23:425–434.
  * 28 Williams CJ, Headd JJ, Moriarty NW, Prisant MG, Videau LL, Deis LN, Verma V, Keedy DA, Hintze BJ, Chen VB, Jain S, Lewis SM, Arendall WB, Snoeyink J, Adams PD, Lovell SC, Richardson JS, Richardson DC (2018) Molprobity: More and better reference data for improved all-atom structure validation. Protein Sci 27:293–315.
  * 29 Spronk CAEM, Nabuurs SB, Krieger E, Vriend G, Vuister GW (2004) Validation of protein structures derived by NMR spectroscopy. Prog Nucl Magn Reson Spectrosc 45:315–337.
  * 30 Schwieters CD, Kuszewski JJ, Tjandra N, Clore GM (2003) The Xplor-NIH NMR molecular structure determination package. J Magn Reson 160: 65–73.
  * 31 Schwieters CD, Kuszewski JJ, Clore GM (2006) Using Xplor-NIH for NMR molecular structure determination. Prog Nucl Magn Reson Spectrosc 48:47–62.
  * 32 Güntert P, Automated NMR Structure Calculation With CYANA. In: Protein NMR Techniques (2004) Methods in Molecular Biology, vol 278. Humana Press, pp 353-378
  * 33 Fowler NJ, Sljoka A, Williamson MP (2021) The accuracy of NMR protein structures in the Protein Data Bank. Structure 29:1430–1439.
  * 34 Vögeli B (2014) The Nuclear Overhauser Effect from a quantitative perspective. Prog Nucl Magn Reson Spectrosc 78:1–46.
  * 35 Chalmers G, Glushka JN, Foley BL, Woods RJ, Prestegard JH (2016) Direct NOE simulation from long MD trajectories. J Magn Reson 265:1–9.
  * 36 Gilpin W (2015) PyPDB: a Python API for the Protein Data Bank. Bioinformatics 32:159–160.
  * 37 Steinegger M, Söding J (2017) MMseqs2 enables sensitive protein sequence searching for the analysis of massive data sets. Nat Biotechnol 35:1026–1028.
  * 38 S. J. Hubbard and J. M. Thornton. Naccess, 1993.
  * 39 Lee B, Richards FM (1971) The interpretation of protein structures: Estimation of static accessibility. J Mol Biol 55:379–400.
  * 40 Zhou AQ, O’Hern CS, Regan L (2012) The power of hard-sphere models: Explaining side-chain dihedral angle distributions of Thr and Val. Biophys J 102:2345–2352.
  * 41 Zhou AQ, O’Hern CS, Regan L (2014) Predicting the side-chain dihedral angle distributions of nonpolar, aromatic, and polar amino acids using hard sphere models. Proteins Struct Funct Bioinf 82:2574–2584.
  * 42 Caballero D, Jukka Määttä, Zhou AQ, Sammalkorpi M, O’Hern CS, Regan L (2014) Intrinsic $\alpha$-helical and $\beta$-sheet conformational preferences: A computational case study of alanine. Protein Sci 23:970–980.
  * 43 Weis S, Schonhoferer PWA, Schaller FM, Schroter M, Schroder-Turk GE (2017) Pomelo: a tool for computing Generic Set Voronoi Diagrams of Aspherical Particles of Arbitrary Shape. EPJ Web of Conferences 140:06007
  * 44 Abraham MJ, Murtola T, Schulz R, Páll S, Smith JC, Hess B, Lindahl E (2015) GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers. SoftwareX 1-2:19–25.
  * 45 Jorgensen WL (1981) Transferable intermolecular potential functions for water, alcohols, and ethers. Application to liquid water. J Am Chem Soc 103:335–340.
  * 46 Best RB, Hummer G (2009) Optimized molecular dynamics force fields applied to the helix-coil transition of polypeptides. J Phys Chem B, 113:9004–9015.
  * 47 Lindorff-Larsen K, Piana S, Palmo K, Maragakis P, Klepeis JL, Dror RO, Shaw DE (2010) Improved side-chain torsion potentials for the Amber ff99SB protein force field. Proteins Struct Funct Bioinf 78:1950–1958.