Statistical Physics of the Symmetric Group

TL;DR

This paper explores the statistical physics of systems modeled by permutations of ordered lists, revealing complex phase behavior with multiple regimes and transitions, and discusses extensions to more complex combinatorial state spaces.

Contribution

It introduces a novel analysis of the symmetric group as a state space in statistical physics, uncovering unique phase transitions and regimes.

Findings

Identification of five physical regimes with two transition temperatures

Discovery of a triple point and a quadruple point in the phase diagram

Extension framework for more complex combinatorial state spaces

Abstract

Ordered chains (such as chains of amino acids) are ubiquitous in biological cells, and these chains perform specific functions contingent on the sequence of their components. Using the existence and general properties of such sequences as a theoretical motivation, we study the statistical physics of systems whose state space is defined by the possible permutations of an ordered list, i.e., the symmetric group, and whose energy is a function of how certain permutations deviate from some chosen correct ordering. Such a non-factorizable state space is quite different from the state spaces typically considered in statistical physics systems and consequently has novel behavior in systems with interacting and even non-interacting Hamiltonians. Various parameter choices of a mean-field model reveal the system to contain five different physical regimes defined by two transition temperatures, a triple point, and a quadruple point. Finally, we conclude by discussing how the general analysis can be extended to state spaces with more complex combinatorial properties and to other standard questions of statistical mechanics models.

Tables

Table 1. Table 1 : Functions defining boundaries between parameter regimes.

Regime Transition	Boundary in Parameter Space
Order/Partial-Order	${\lambda }_1=\frac{1}{\beta }\ln N$
Order/Order-Partial-Order Metastability	${\lambda }_2=\frac{1}{\beta }{W}_{-1}\left(-\frac{N}{e}{e}^{-\beta {\lambda }_1}\right)$
Order/Disorder	${\lambda }_2=-\frac{{\lambda }_1}{1-1/N}$

Equations

$$Z=\sum _{k=1}^{N!}{e}^{-\beta {ϵ}_k},$$

$$\vec{\omega }=({\omega }_1,\dots ,{\omega }_N).$$

$$Sym(\omega ):=\text{Set of All Permutations of (ω1,…,ωN)}.$$

$${\mathcal{H}}_N(\{{\theta }_i\})=\sum _{i=1}^N{\lambda }_i{I}_{{\theta }_i\ne {\omega }_i},$$

$$I_A\equiv \{\begin{array}{ll}1& \text{if A is true}\\ 0& \text{if A is false}\end{array}.$$

$$Z_N(\{\beta {\lambda }_i\})=\sum _{\vec{\theta }\in Sym(\omega )}\text{exp}\left(-\beta \sum _{i=1}^N{\lambda }_i{I}_{{\theta }_i\ne {\omega }_i}\right),$$

$$Z_N(\{\beta {\lambda }_i\})=\sum _{j=0}^N{d}_j{\Pi }_j\left(e^{-\beta {\lambda }_1},\dots ,e^{-\beta {\lambda }_N}\right),$$

$${\Pi }_k(x_1,\dots ,x_N)=\frac{1}{k!}{\left[\frac{d^k}{d{q}^k}\prod _{i=1}^N(1+q{x}_i)\right]}_{q=0}.$$

$$Z_N(\{\beta {\lambda }_i\})$$

$$\displaystyle=\int^{\infty}_{0}ds\,e^{-s}\prod_{\ell=1}^{N}\Big{[}1+(s-1)e^{-\beta\lambda_{\ell}}\Big{]},$$

$$Z_N(\beta {\lambda }_0)$$

$$={\int }_0^{\infty }ds{e}^{-s}{\left(1+(s-1)e^{-\beta {\lambda }_0}\right)}^N$$

$$j\equiv \sum _{i=1}^N{I}_{{\theta }_i\ne {\omega }_i},$$

$$g_N(j)=\left(\genfrac{}{}{0ex}{}{N}{j}\right)d_j$$

$$\beta F$$

$$\simeq N\beta {\lambda }_0-\left(e^{\beta {\lambda }_0}-N-1\right)+F_0(N),$$

$$⟨j⟩\simeq N-e^{\beta {\lambda }_0}.$$

$$k_B{T}_c=\frac{{\lambda }_0}{\ln N}.$$

$${\mathcal{H}}_N(\{{\theta }_i\})={\lambda }_0\sum _{i=1}^N{I}_{{\theta }_i\ne {\omega }_i}.$$

$${\mathcal{H}}_{\text{system}}={\mathcal{H}}_1\oplus {\mathcal{H}}_2\oplus \cdots \oplus {\mathcal{H}}_N,$$

$${\mathcal{S}}_{\text{system}}\ne {\mathcal{S}}_1\otimes {\mathcal{S}}_2\otimes \dots \otimes {\mathcal{S}}_N.$$

$$\Delta F(j)$$

$$\mathcal{H}(\{{\theta }_i\})=\sum _i{\lambda }_i{I}_{{\theta }_i\ne {\omega }_i}.$$

$$\mathcal{H}(\{{\theta }_i\})=\sum _i{\lambda }_i{I}_{{\theta }_i\ne {\omega }_i}+\frac{1}{2}\sum _{i,j}{\mu }_{ij}{I}_{{\theta }_i\ne {\omega }_i}{I}_{{\theta }_i\ne {\omega }_i}+\cdots .$$

$$Z_N(\beta ;{\lambda }_1,{\lambda }_2)$$

$$\frac{\beta {\lambda }_2}{2N}\sum _{i,j=1}^N{I}_{{\theta }_i\ne {\omega }_i}{I}_{{\theta }_j\ne {\omega }_j}),$$

$$Z_N(\beta ;{\lambda }_1,{\lambda }_2)=\sum _{j=0}^N{g}_N(j)e^{-\beta \mathcal{E}(j)},$$

$$\mathcal{E}(j)={\lambda }_{1}j+\frac{{\lambda }_2}{2N}{j}^2$$

$$F_N(j,\beta )={\lambda }_{1}j+\frac{{\lambda }_2}{2N}{j}^2-\frac{1}{\beta }\ln g_N(j).$$

$$g_N(j)=\frac{\Gamma (N+1)}{\Gamma (j+1)\Gamma (N-j+1)}\frac{\Gamma (j+1,-1)}{e},$$

$$f_N(j,\beta )={\lambda }_{1}j+\frac{{\lambda }_2}{2N}{j}^2+\frac{1}{\beta }\ln \Gamma (N-j+1)+f_0$$

$$\frac{\partial }{\partial j}{f}_N(j,\beta )={\lambda }_1+\frac{{\lambda }_2}{N}j-\frac{1}{\beta }{\psi }_0(N-j+1)=0.$$

Full text

I Introduction

Chains of amino acids are important components of biological cells, and for such chains the specific ordering of the amino acids is often so fundamental to the resulting function and stability of the folded chain that if major deviations from the correct ordering were to occur, the final chain could fail to perform its requisite function within the cell, proving fatal to the organism.

More specifically, we see the relevance of correct ordering in the study of protein structure, which is often divided into the protein folding and protein design problem. While the protein folding problem concerns finding the three-dimensional structure associated with a given amino acid sequence, the protein design problem (also termed the inverse-folding problem; see Figure 1) concerns finding the correct amino acid sequence associated with a given protein structure.

An aspect of one solution to the protein design problem is to maximize the energy difference between the low-energy folded native structure and the higher energy misfolded/denatured structures. In doing so, one takes native structure as fixed and then determines the sequence yielding the minimum energy, under the assumption (termed the ”fixed amino-acid composition” assumption) that only certain quantities of amino-acids appear in the chain Morrissey and Shakhnovich (1996). In this resolution (specifically termed heteropolymer models Shakhnovich and Gutin (1993a) Shakhnovich and Gutin (1993b)) the correct amino acid sequence is found by implementing an MC algorithm in sequence space given a certain fixed amino acid composition. This entails assuming the number of various types of amino acids does not change, and distinct states in sequence space are permutations of one another. For example, for a polypeptide chain with $N$ residues, rather than searching over the entire sequence space (of size $20^{N}$), one searches over a space of sequences (of size $N!/n_{1}!n_{2}!\ldots n_{20}!$) which are defined by a fixed number of each amino acid.

This aspect of the protein design problem alerts one to a gap in the statistical mechanics literature. Namely, there do not seem to be any simple and analytically soluble statistical mechanics models where the space of states is defined by permutations of a list of components.

We can take steps toward constructing such a model by considering reasonable general properties it should have. If we assume there was a specific sequence of components which defined the lowest energy sequence and was thermodynamically stable in the model, then deviations from this sequence would be less stable. Because of the role sequences of molecules play in biological systems, it is worth asking what features we expect such sequences to have from the perspective of modeling in statistical mechanics.

In Section II we introduce the model, and compute an exact partition function which displays what we term “quasi”-phase transitions—a transition in which the sequence of lowest energy becomes entropically disfavored above a certain temperature. In Section III, we extend the previous model by adding a quadratic mean field interaction term and show that the resulting system displays two transition temperatures, a triple point, and a quadruple point. In Section IV, we discuss various ways we can extend this model in theoretical or more phenomenological directions.

II System and Partition Function

Our larger goal is to study equilibrium thermodynamics for a system defined by permutations of a set of $N$ components where each unique permutation is defined by a specific energy. In general, we should consider the case where the set of $N$ components consists of $L$ types of components for which if $n_{k}$ is the number of repeated components of type $k$, then $\sum_{k=1}^{L}n_{k}=N$. For simplicity, however, we will take $n_{k}=1$ for all $k$ so that each component is of a unique type and $L=N$.

To study the equilibrium thermodynamics of such a system with a fixed $N$ at a fixed temperature $T$, we need to compute its partition function. For example, for a sequence with $N$ components (with no components repeated), there are $N!$ microstates the system can occupy and assuming we label each state $k=1,\ldots,N!-1,N!$, and associate an energy $\epsilon_{k}$ with each state, then the partition function would be

[TABLE]

where $\epsilon_{k}$ for each state $k$ could be reasoned from a more precise microscopic theory of how the components interact with one another. Phenomenologically, Eq.(1) would be the most precise way to construct a model to study the equilibrium properties of permutations, but because it bears no clear mathematical structure, it is unenlightening from a theoretical perspective.

Instead we will postulate a less precise, but theoretically more interesting model. For most ordered chains in biological cells, there is a single sequence of components which is the “correct” sequence for a particular macrostructure. Deviations from this correct sequence are often disfavored because they form less stable macrostructures or they fail to perform the original function of the “correct” sequence. With the general properties of such sequences in mind, we will abstractly represent our system as consisting of $N$ sites which are filled with particular coordinate values denoted by $\omega_{k}$. That is, we have an arbitrary but fixed coordinate vector $\vec{\omega}$ expressed in component form as

[TABLE]

We will take the collection of components $\{\omega_{k}\}$ as intrinsic to our system, and thus take the state space of our system to be the set of all the vectors whose ordering of components can be obtained by permuting the components of $\vec{\omega}$, i.e., all permutations of $\omega_{1},\ldots,\omega_{N}$. We represent an arbitrary state in this state space as $\vec{\theta}=(\theta_{1},\ldots,\theta_{N})$, where the $\theta_{k}$ are drawn without repeat from $\{\omega_{k}\}$. Formally, we would say our space of states is isomorphic to the symmetric group on $\vec{\omega}$ (Dixon and Mortimer (1996)). We will thus denote our state space as

[TABLE]

and then an arbitrary state $\vec{\theta}$ is just an element element of this set.

As a first formulation of the model, we will take $\vec{\theta}_{0}=\vec{\omega}$ (the correct sequence) to represent the zero energy state in the system, and for each component $\theta_{i}$ of an arbitrary vector $\vec{\theta}$ which differs from the corresponding component $\omega_{i}$ in $\vec{\omega}$, there is an energy cost of $\lambda_{i}>0$. The Hamiltonian is then

[TABLE]

where $\theta_{i}$ and $\omega_{i}$ are components of vectors $\vec{\theta}$ and $\vec{\omega}$ respectively, and $I$ is defined by

[TABLE]

We note that although we label our general state as $\vec{\theta}=(\theta_{1},\ldots,\theta_{N})$, the components $\theta_{1},\ldots,\theta_{N}$ can only take on mutually-exclusive values from the set $\{\omega_{k}\}$.

We want to explore the equilibrium thermodynamics of a system with a Hamiltonian of Eq.(4). This amounts to calculating the partition function

[TABLE]

where $Sym(\omega)$ is again the set of all permutations of the components of $(\omega_{1},\ldots,\omega_{N})$. To find a closed form expression for the partition function, we group terms in Eq.(6) according to the number of ways to completely reorder $j$ components in $\vec{\omega}$ while keeping the remaining components fixed. Each such reordering (i.e., each value of $j$) is associated with a sum over products of $e^{-\beta\lambda_{i}}$ terms with $j$ factors of $e^{-\beta\lambda_{i}}$ (for various $i$) in each term. The total partition function is a sum of all such reorderings for all $j$s from [math] to $N$. As can be seen from a direct expansion of Eq.(6), we have

[TABLE]

where $d_{j}$, termed the number of derangements of a list of $j$ (Chuan-Chong and Khee-Meng (1992)), is the number of ways to completely reorder a list of $j$ elements. The quantity $\Pi_{j}(x_{1},\ldots,x_{N})$, termed the $j$th elementary symmetric polynomial on $n$ (Borwein and Erdélyi (2012)), is the sum of all ways to multiply $j$ elements out of the $N$ term set $\{x_{1},\ldots,x_{N}\}$. For example, $\Pi_{2}(x_{1},x_{2},x_{3})=x_{1}x_{2}+x_{2}x_{3}+x_{3}x_{1}$. By definition

[TABLE]

By the definition of the incomplete gamma function as $\Gamma(x,\alpha)=\int^{\infty}_{\alpha}dt\,t^{x-1}e^{-t}$ and its relation to derangements (i.e., $d_{j}=\Gamma(j+1,-1)/e$, see Weisstein (2002a)), we then find

[TABLE]

which is the desired closed-form expression for the partition function of this system.

With Eq.(9), the problem of abstractly studying a thermal system of permutations with Hamiltonian Eq.(4) is, from the perspective of equilibrium statistical mechanics, now complete. However, there are still some physical and theoretical results which can be teased from this formalism. Specifically, we can ask whether this system exhibits phase transitions. To answer this question, it would prove more analytically tractable to take $\lambda_{i}=\lambda_{0}$ for all $i$. With this condition, our partition function simplifies to

[TABLE]

where we transformed our Hamiltonian as ${\cal H}_{N}(\{\theta_{i}\})\to{\cal H}(j)=\lambda_{0}j$ with $j$, defined as

[TABLE]

the number of components of $\vec{\theta}$ which are not equal to the corresponding component in $\vec{\omega}$. We call $j$ the number of incorrect components of $\vec{\theta}$, and if $j=N$ we say $\vec{\theta}$ is completely disordered. The factor $g_{N}(j)$ in Eq.(10) is defined as

[TABLE]

is the number of ways to reorder a list of $N$ elements so that $j$ elements are no longer in their original position. This combinatorial definition of $g_{N}(j)$ will prove useful when we explore the phase behavior of more complex models of permutations.

From the form of Eq.(10), it is clear that, physically, its associated Hamiltonian is not realistic as it places distinct permutations (which in any true physical system most likely have quite different energy properties) in the same degenerate energy state. Still, from a theoretical perspective, the simplicity of this model makes it a suitable starting point for studying the general properties of systems of permutations.

II.1 Phase-like Behavior of “Non-Interacting” System

We can investigate the phase-like behavior of the system defined by the Hamiltonian Eq.(4) (for constant $\lambda_{i}$ across $i$), by applying steepest descent to Eq.(11) in the $N\gg 1$ limit. Doing so, we find we can approximate the free energy of the system to be

[TABLE]

where $F_{0}(N)$ is a $\beta\lambda_{0}$ independent constant. Noting that $\langle j\rangle=-\partial\ln Z_{N}(\beta\lambda_{0})/\partial(\beta\lambda_{0})=\partial F/\partial\lambda_{0}$, we find the average number of incorrect components satisfies the following equation of state:

[TABLE]

By Eq.(12), we can infer that $\langle j\rangle$ must be greater than or equal to 0. However, the right-hand-side of Eq.(15) exhibits no such explicit constraint. Thus we can infer there is a phase-like transition in our system at the temperature

[TABLE]

Below this temperature, we must have $\langle j\rangle\simeq 0$ and thus the “correct permutation” has the lowest free energy and is thermodynamically favored; above this temperature, $\langle j\rangle>0$ and the system is in a disordered phase where the previous lowest energy “correct-permutation” is energetically disfavored.

Interestingly, this transition arises from the naively non-interacting Hamiltonian

[TABLE]

We say “naively non-interacting” because Eq.(17) consists of a sum over linear functions of a single index $i$, and thus does not suggest any coupling between terms of differing index. However, statistical mechanics tells us that the energy of a system isn’t the only thing which determines the thermodynamic behavior of a system. Indeed we have to consider entropic contributions as well, and in this system the entropy (as it is a function of $j$) can drive thermodynamic behavior. In other words, although the Hamiltonian is depicted as non-interacting and can set-theoretically be represented as

[TABLE]

our system really exhibits interactions between components because our total space of states ${\cal S}$ cannot be factorized:

[TABLE]

Thus the “non-interacting” system exhibits a transition at Eq.(16) due to the coupled nature of the state space. As discussed in the subsequent section, we term this transition a “quasi-phase transition” because it does not bear all of the standard properties we expect in phase transitions.

II.1.1 Not a True Phase Transition

We claim the system does not exhibit true phase transition behavior because many of these results are not consistent with the traditional thermodynamic definition of phase transitions. For one, phase transitions are associated with divergences in the derivatives of the free energy, but there is no divergence in the free energy associated with the partition function Eq.(11) for possible parameter values. Also, the result Eq.(15) arises from the steepest descent approximation which makes $\langle j\rangle$’s temperature dependence near $\langle j\rangle=0$ appear non-differentiable when by Eq.(11) it is actually differentiable over its entire domain. Finally, with Eq.(10) we can define a Landau free energy $F(j)$ for this system according to $Z=\sum_{j}e^{-\beta F(j)}$, and what we may ordinarily label as a phase transition (i.e., going from $\langle j\rangle=0$ to $\langle j\rangle\neq 0$) arises, not from changes in the functional form of the Landau free energy as we see in real phase transitions, but from changes in the excluded region of the Landau free energy (See Figure 2). Because the functional form of the free energy remains the same we observe no true phase transition.

Alternatively, a heuristic argument for the non-existence of phase transitions in our permutation model is mathematically very similar to the Landau argument (Landau et al. (1980)) for the non-existence of transitions in 1d Ising Models. For our permutation system with $N$ lattice sites, the state of zero energy and zero entropy consists of every site being occupied by its correct component. To increase the energy of this system, we can choose $j$ sites to contain incorrect components, thus giving us an energy ${\cal H}_{j}=\lambda_{0}j$. The number of ways we can choose these $j$ components is given by Eq.(13) Thus, upon introducing $j\neq 0$ incorrect components, the change in the Landau free energy of our system is

[TABLE]

where we took these results in the $1\ll j\ll N$ limit and used $d_{j}\simeq j!/e$. In the thermodynamic ($N\to\infty$) limit, we find that $\Delta F(j)\to-\infty$ meaning there is no non-zero $T$ at which the entropic contribution becomes subdominant to the energy. Thus the system exhibits no phase transition.

III Partition Function for Interacting Model

When we first considered a model of thermal permutations, we began with a Hamiltonian where sites did not interact with one another and each had a site-dependent energy cost for being incorrectly occupied:

[TABLE]

More generally, we can consider Hamiltonians with an arbitrary number of multiple-site interaction terms. Such a Hamiltonian could be written as

[TABLE]

The first term in Eq.(22) associates an energy cost of $\lambda_{i}$ with incorrectly occupying the component at position $i$. The second term models interactions between sites where the correct (or incorrect) occupation of a single site determines the energy of another. The exact values of $\lambda_{i}$ and $\mu_{ij}$ could be chosen to ensure the “correct” state ($\theta_{i}=\omega_{i}$ for all $i$) is non-degenerate as in the non-interacting model. The ellipses represent higher order interactions in this framework. Hamiltonians such as Eq.(22) should be more physically relevant as they would correspond to systems where the energy cost for deviating from the lowest energy permutation is not simply linear but could be represented as a tensor valued fitting function.

We can make progress in studying the thermodynamics of more general Hamiltonians like Eq.(22) by first only considering first- and second-order interaction terms and taking the interactions to be constants: $\lambda_{i}=\lambda_{1}$ for all $i$; $\mu_{ij}=\lambda_{2}/N$ for all $i,j$. The factor of $1/N$ is chosen so that the second term matches the extensive scaling of the first term. The partition function for such parameter selections is then

[TABLE]

where $\lambda_{1}$ and $\lambda_{2}$ are interaction parameters with units of energy. We can also write this partition function in the Eq.(12) basis as

[TABLE]

where $g_{N}(j)$ is defined in Eq.(13) and

[TABLE]

is the energy function for the system.

III.1 Calculating Order Parameter

Our goal is to analyze the “quasi”-phase behavior of this system in a way analogous to our analysis for the non-interacting system. To do so we begin with the Landau free energy function

[TABLE]

Alternative starting points for this derivation are presented in Appendix A. Our system is constitutively discrete, so it is not precisely correct to discuss our free energy in the language of analysis, but given our expression for Eq.(13) we can map this system to a continuous one which bears the same thermodynamic properties and for which analysis is appropriate. Specifically, if we take $j$ to be continuous and use the identity $\Gamma(x+1)=x!$, we can write

[TABLE]

With the approximation $\Gamma(j+1,-1)\simeq\Gamma(j+1)$ and the substitution Eq.(27), Eq.(26) then becomes

[TABLE]

where we defined our approximated free energy as $f_{N}(j,\beta)$ and collected the $j$ independent constants into $f_{0}$. Now Eq.(28) is fully continuous and amenable to analysis. To find the thermodynamic equilibrium of this system, we need to find the value of $j$ for which $\partial f_{N}(j,\beta)/\partial j=0$ and $\partial^{2}f_{N}(j,\beta)/\partial j^{2}>0$. For the first condition we have

[TABLE]

For $x\geq 0.6$ we have

[TABLE]

as can be affirmed by Taylor expansion or plots of each side. Since the argument of $\psi_{0}(N-j+1)$ is bounded below by 1, the approximation in Eq.(30) can be applied to Eq.(29). With this substitution, and setting the result to be valid for the equilibrium value $j=\overline{j}$, we then find the constraint

[TABLE]

which has the solution

[TABLE]

where $W$ is the (branch unspecified) Lambert W function, defined by

[TABLE]

To specify the branch of the $W$ which corresponds to a stable equilibrium we compute the second derivative of our free energy at this derived critical point. Doing so yields

[TABLE]

Thus Eq.(32) (for $\lambda_{2}>0$) yields a free energy minimum for

[TABLE]

and yields a maximum for the inverse condition. This amounts to stating that the stable equilibrium for $\overline{j}$ is defined by the principal branch of the Lambert W function where $W=W_{0}\geq-1$, and the unstable equilibrium for $\overline{j}$ is defined by the negative branch where $W=W_{-1}<-1$.

Thus, the order parameter for this system is

[TABLE]

We note that taking $\lambda_{2}\to 0$ and using $W(x)=x+{\cal O}(x^{2})$ for $|x|\ll 1$ returns us to the non-interacting result Eq.(15).

For completeness, we define the value of $j$ which yields a free energy maximum as $\overline{j}_{-1}$; it is related to Eq.(36) by replacing the principal branch function $W_{0}$ with $W_{-1}$.

III.2 Discussion of Parameter Space

In the previous section, we found that the order parameter for this system was given by Eq.(36). We noted that this solution represented a local minimum of the free energy as long as the Lambert W function satisfied $W=W_{0}>-1$. Thus when this condition is violated, $\overline{j}_{0}$ is no longer a valid stable equilibrium, and our system has undergone a “quasi”-phase transition or simply a transition.

Moreover, our values of $j$ are bounded below by $j=0$ and bounded above by $j=N$, neither conditions of which are naturally constrained by Eq.(36). Thus these two conditions are associated with two other transitions. In all, then, there are three conditions which define the quasi-phase boundaries in this system.

While there are three conditions which define transitions in this system, there are in fact five distinct regimes of parameter space. We can obtain a qualitative sense of these regimes by creating schematic plots of the free energy Eq.(28) for various parameter values of $\lambda_{1}$ and $\lambda_{2}$. The possible plots can be placed into five categories according to the plot’s stable or metastable $j$ values. We depict these possible plots in Figure 3. We note that only the free energy plots with valid values of $\overline{j}_{0}$ contain what we normally consider a thermodynamic equilibrium; the other plots have “stable” values of $j$ arising only from the $j=0$ and/or $j=N$ boundary conditions.

Qualitatively, we can name the states according to the sequence space to which their equilibrium values of $j$ correspond. We know for $j=0$, our system is in a state with zero incorrect components in $\vec{\theta}$ and hence the system is “perfectly ordered” or just ”ordered”. Conversely for $j=N$ our system has $N$ incorrect components and hence the system is “completely disordered” or just “disordered”. The in-between case of $j=\overline{j}$ where $0<\overline{j}<N$ can be given the related label of “partially-ordered”. Thus, the regime names associated with our possible values of the order parameters are

• •

Ordered Regime ($\boldsymbol{j=0}$ stable): Neither $\overline{j}_{0}$ or $\overline{j}_{-1}$ exist; $f(N,\beta)>0$.

• •

Disordered Regime **($\boldsymbol{j=N}$ stable) **: Neither $\overline{j}_{0}$ or $\overline{j}_{-1}$ exist; $f(N,\beta)<0$.

• •

Partially-Ordered Regime ($\boldsymbol{j=\overline{j}_{0}}$ stable): Only $\overline{j}_{0}$ exists.

• •

Order and Disorder Metastable Regime ($\boldsymbol{j=0}$ and $\boldsymbol{j=N}$ stable): Only $\overline{j}_{-1}$ exists.

• •

Order and Partial-Order Metastable Regime **($\boldsymbol{j=0}$ and $\boldsymbol{j=\overline{j}_{0}}$ stable) **: Both $\overline{j}_{0}$ and $\overline{j}_{-1}$ exist.

We note that it seems to be a fundamental feature (or a lack of one) of this system, that the free energy Eq.(28) does not admit a metastability between partial-order and disorder.

III.3 Monte-Carlo Generated Parameter Space

With these regime definitions, we can depict the parameter space graphically. In Figure 3 we showed the possible forms of the free energy for this system where each was categorized according to the existence of the local minimum critical point $\overline{j}_{0}$, the existence of the local maximum critical point $\overline{j}_{-1}$, and the sign of the quantity $f_{N}(j,\beta)$. We can extrapolate this categorization to $\lambda_{1}-\lambda_{2}$ parameter space, by determining which regions of parameter space correspond to specific plots in Eq.(28). Doing so through the Monte Carlo procedure described in Appendix B, we generated 10,000 points of the parameter space diagram in Figure 4 for $\beta$ set to 1. We note that the parameter space exhibits five regimes separated by three lines cited in Table 1 each of which correspond to the three conditions mentioned at the beginning of this section. These lines can be derived analytically (as shown in Appendix C) by considering the conditions in turn and which regimes they serve to connect.

III.4 Triple and Quadruple Points and Transition Temperatures

From Figure 4 we see that our system is characterized by two points where there is a coexistence between multiple regimes. Given that $W_{-1}(-e^{-1})=-1$, we have that the Ordered, Partially-Ordered, and Order-Partially-Ordered Metastability coexistence point is characterized by the condition

[TABLE]

These conditions characterize the system’s triple point.

Similarly, for $N\gg 1$, the Partially-Ordered, Disordered, Order-Partially-Ordered Metastability, and Order-Disorder Metastability coexistence point is characterized by the condition

[TABLE]

This condition characterizes the quadruple point of the system.

Figure 4 also depicts the possible regimes of our system for a given temperature and various Hamiltonian parameters $\lambda_{1}$ and $\lambda_{2}$. More physically, we may be interested in knowing what are the “quasi” phase properties of a system with a fixed $\lambda_{1}$ and $\lambda_{2}$ and a variable temperature. That is, what are the temperatures which define the various transitions between regimes in the system?

An arbitrary permutation system for a fixed $N$ and at a variable temperature is characterized by a specific energy function Eq.(25). Such a system is therefore defined by a specific $\lambda_{1}$ and $\lambda_{2}$, and the system can be associated with a particular point (and hence region) in the parameter space of Figure 4. As we vary the temperature of this system, the temperature dependent regime-coexistence lines change and if they change in such a way as to extend the region of a regime to newly encompass our original point then our system has undergone a transition. In this way, we can define the temperatures which characterize various possible transitions of this system.

First, from Figure 4 and Eq.(85) we note the regime-coexistence line between the partially-ordered and disordered regime is independent of temperature, and so there is no critical temperature defining a partial-order to disorder transition.

From Eq.(78), we can infer that the partial-order to order transition is characterized by moving below the temperature

[TABLE]

Contingent on which region of parameter space the system lies, this temperature also characterizes the disorder to order-disorder metastability transition and the partial-order to order-partial-order metastability transition.

And from Eq.(81), we can solve for the associated transition temperature given fixed $\lambda_{1}$ and $\lambda_{2}$ to find

[TABLE]

where this expression is only relevant for $-\lambda_{1}<\lambda_{2}<0$ and $\lambda_{1}>k_{B}T\ln N$. Moving above this temperature leads to the order to order-partial-order metastability transition.

IV Discussion

In this work, motivated by an abstraction of a foundational problem in protein design, we posited and analyzed the basic properties of a statistical physics model of permutations. Formally, we considered a simple statistical physics model where the space of states for $N$ lattice sites was isomorphic to the symmetric group of degree $N$ Dixon and Mortimer (1996), and where the energy of each permutation was a function of how much the permutation deviates from the identity permutation.

In this model, we found that due to a state space which could not be factorized in a basis defined by lattice sites, even the superficially non-interacting system can exhibit phase-like transitions, i.e., temperature dependent changes in the value of the order parameter which do note exhibit the properties typically associated which phase transitions in infinite systems. When interactions are introduced through a quadratic mean field term, the system is capable of exhibiting five regimes of thermal behavior, and is characterized by two transition-temperatures corresponding to various quasi-phase transitions.

The introduced model provides us with a basic exactly soluble system for certain interaction assumptions and thus provides a concrete model-based understanding of a system with a non-factorizable state space. Because of its utility and the type of results obtained, the model deserves to be subject to the standard extensions of typical canonical models in statistical mechanics. In particular we hope to extend it to non-trivial site dependent interactions. For example, a nearest neighbor interaction Hamiltonian of the kind which characterize the Ising Model,

[TABLE]

would be an alternative physical extreme to the mean-field interactions considered in Section III.

We could also consider a generalized chain of components where the interactions between sites or the cost for an incorrectly filled site is not constant but is drawn from a distribution of values. Such a system of quenched disorder would characterize a permutation glass which may contain interesting results due to the unique nature of the state space.

Supposing it is possible to define more interesting interactions models, a natural investigation would concern the renormalization group properties of the system. Specifically, we would be interested in how would one sum over specific states (as characteristic of a renormalization group transformation) when the state space of a system looks like,

[TABLE]

i.e., is not factorizable along lattice sites.

Finally, to connect this model of permutations to problems more relevant to protein design it would prove necessary to incorporate the possibility of repeated components or the background geometry of a lattice chain.

Acknowledgements.

MW thanks Verena Kaynig-Fittkau for her unpublished work that inspired this investigation and Amy Gilson, Abigail Plummer, Vinothan Manoharan, and Michael Brenner for helpful discussions.

Appendix A Alternative Derivations of Eq.(31)

A.1 Hubbard-Stratonovich Approach

We re-derive Eq.(31) using the Hubbard-Stratonovich method. We start this derivation assuming $\lambda_{2}=-|\lambda_{2}|$; we will later see our resulting free energy can be analytically continued to the $\lambda_{2}>0$ case.

For $\lambda_{2}=-|\lambda_{2}|$, the partition function is

[TABLE]

Then, applying the identity

[TABLE]

we have

[TABLE]

where $Z_{N}(x)\equiv\sum_{j=0}^{N}g_{N}(j)e^{-jx}$. From Eq.(11) we found

[TABLE]

So Eq.(45) becomes

[TABLE]

The function to which we apply steepest descent is then

[TABLE]

Computing the conditions for $\partial_{s}h(s=\overline{s},x=\overline{x})=0$ and $\partial_{x}h(s=\overline{s},x=\overline{x})=0$ we obtain, respectively,

[TABLE]

Solving for $\overline{s}$ in the first equation we have

[TABLE]

and with the second equation we have the condition

[TABLE]

Substituting the second condition into the first yields

[TABLE]

or

[TABLE]

Thus, the solution for $\overline{s}$ can be expressed in terms of the Lambert W function as

[TABLE]

For $\lambda_{2}>0$, we can employ the complex version of the Hubbard-Stratonovich identity:

[TABLE]

Working through an analogous steepest descent procedure, we find that the equilibrium value for $\overline{s}$ is

[TABLE]

which could have been extrapolated from Eq.(55) by taking $|\lambda_{2}|\to-\lambda_{2}$. From this expression for the equilibrium condition, and as an analogy with the non-interacting case, it turns out the order parameter in this case is not $\overline{s}$ but rather $\overline{s}-1$.

A.2 Gibbs-Bogoliubov Inequality Derivation of Eq.(31)

We re-derive Eq.(31) using the Gibbs-Bogoliubov InequalityYeomans (1992). The inequality is

[TABLE]

The Hamiltonian which defines our system is

[TABLE]

and our variational Hamiltonian is instead

[TABLE]

From Eq.(11) we know

[TABLE]

We can also define

[TABLE]

as the average with respect to our variational Hamiltonian Eq.(60). Thus, Eq.(58) becomes

[TABLE]

Differentiating $f$ with respect to $\lambda_{0}$ allows us to compute the maximum of this quantity. Given $\langle j\rangle_{0,\,N}=-\partial\ln Z_{N}(\beta\lambda_{0})/\partial(\beta\lambda_{0})$, we then find

[TABLE]

which if we take to be zero at some $\lambda_{0}={\overline{\lambda}_{0}}$ gives us

[TABLE]

To compute these derivatives we make use of various identities. First we note

[TABLE]

so

[TABLE]

This last equality implies

[TABLE]

and so Eq.(65) becomes

[TABLE]

To compute these quantities we need to approximate the partition function for our variational system. Using the method of steepest descent

[TABLE]

and so we have

[TABLE]

Computing the relevant derivatives in Eq.(69) we find

[TABLE]

Neglecting subleading terms of ${\cal O}(1/N)$ (a choice only valid for $\langle j\rangle_{0,N}\gg 1$), solving for ${\overline{\lambda}_{0}}$, and using Eq.(71) we then obtain the equilibrium constraint

[TABLE]

which when solved for $\langle j\rangle_{0,N}/N$ gives us

[TABLE]

or, given our approximations and limiting expressions, the result

[TABLE]

Appendix B Monte-Carlo Procedure for Parameter Space

To generate Figure 4, we implemented the following MC algorithm:

Uniformly sample two points for $\lambda_{1}$ and $\lambda_{2}$ separately from within a certain bounded domain. 2. 2.

Draw the free energy curve Eq.(28) corresponding to the sampled values $(\lambda_{1},\lambda_{2})$. 3. 3.

Label the curve according to which schematic curve in Figure 3 it corresponds (i.e., according to its $\overline{j}_{0}$, $\overline{j}_{-}$ and $f_{N}(N,\beta)$ properties). 4. 4.

Color the point to signify the label.

We repeated this procedure for $10,000$ points with $\beta=1$. The regime separation lines were included after the MC procedure from the analytic forms cited in the text.

Appendix C Analytic Functions of Regime Boundaries

C.0.1 Order to Partial-Order Transition

The regime boundary which separates the ordered and the partially ordered regime is defined by the condition $\overline{j}_{0}\geq 0$. For this regime boundary we have the condition

[TABLE]

or

[TABLE]

From a plot of $W_{0}(a\,xe^{x})/x$ for real $a$, we see that $W(axe^{x})/x>1$ if $a>1$ and $W(axe^{x})/x<1$ if $a<1$. Thus this order to partial-order transition is defined by the condition $e^{\beta\lambda_{1}}/N=1$, or

[TABLE]

C.0.2 Order to Order-Partial-Order Metastability Transition

The regime boundary which separates the ordered regime from the order and partial-order metastability regime is defined by the condition $-1\leq W<0$. This condition is where the $\overline{j}_{0}$ and $\overline{j}_{-}$ begin coexisting Weisstein (2002b), thus creating the mutual existence of a local maxima and local minima in Figure 3. Thus for this regime boundary we have the condition

[TABLE]

This condition is valid so long as the argument of $W$ satisfies

[TABLE]

This inequality can only possibly be satisfied for $\lambda_{2}<0$ and if $\lambda_{2}<0$ the right inequality is automatically true. So, our transition condition is given by

[TABLE]

In the Order phase we automatically have $\beta\lambda_{1}>\ln N$, so the LHS of Eq.(81) is greater than or equal to $-e^{-1}$. Moreover, since $\lambda_{2}$ is exclusively negative, at $\beta\lambda_{2}=\ln N$, $\beta\lambda_{1}$ is at a maximum value of $\beta\lambda_{2}=-1$. For $\beta\lambda_{2}\leq-1$, the solution to Eq.(81) is then

[TABLE]

C.0.3 Order to Disorder Transition

The regime boundary which separates the partially ordered regime from the disordered regime is defined by the condition $\overline{j}_{0}\leq N-1$. We set the maximum value of $\overline{j}_{0}$ to $N-1$ rather than $N$ because Eq.(36) is associated with a free energy which diverges at $\overline{j}_{0}=N$ and this approximate result is thus only physical up to $N-1$. Alternatively we could see the maximum condition $\overline{j}_{0}=N-1$ as respecting the fact that Eq.(36) is only valid up to ${\cal O}(N^{-1})$. For this regime boundary we have the condition

[TABLE]

or

[TABLE]

Again, using the condition that $W_{0}(a\,xe^{x})/x>1$ if $a>1$, we find that the critical condition for this transition is $\beta\lambda_{1}+\beta\lambda_{2}-\beta\lambda_{2}/N=0$ or

[TABLE]