text:
# A Solvable Model of Secondary Structure Formation in Random Hetero-Polymers

## 1 Introduction

Proteins are polymeric chains of amino-acids. The successful functioning of a protein in a living organism depends crucially, among other factors, on its ability to fold into a desired three-dimensional structure (its ‘native state’), and to subsequently attach in a very specific way to other macro-molecules. From a biological and medical point of view, it is therefore highly desirable to know which native state corresponds to a given amino-acid sequence, and (conversely, for therapeutic purposes) to know which amino-acid sequence would fold into a desired native state; this requires a quantitative understanding of the physical forces underlying the folding mechanism. A detailed identification of sequence-specific native states will necessarily involve sophisticated (molecular dynamics based) computational methods. However, due to the large number of degrees of freedom of proteins, the complicated nature of the various types of electro-chemical interactions and the so-called ‘hard’ geometric chain constraints of a protein, such computer programmes are unfortunately (as yet) extremely slow. Thus, in order to identify the role and degree of importance of the various folding parameters, a theoretical (i.e. statistical mechanical) analysis would be very welcome.

It is generally assumed that the presently observed population of real proteins has evolved from the larger class of random hetero-polymers, driven by natural selection. This suggests that the study of random hetero-polymers is a natural first step en route towards the statistical mechanical study of proteins. Furthermore, already at an early stage it was recognized , via a theoretical study based on the random energy scheme , that many aspects of protein folding (such as the appearance of ‘mis-folded’ phases, and transitions between folded and unfolded states) can be understood on the basis of equilibrium statistical mechanical calculations for random hetero-polymers. Even simple models with only two types of amino-acids interacting with the water solvent, viz. hydro-phobic amino-acids versus polar ones, can successfully describe the basics of protein folding (see e.g. the so-called HP model ). Further statistical mechanical approaches include replica calculations on polymer chains with Gaussian pair interactions , variational analyses in replica spaces , lattice models  and lattice gas models , to mention but a few. In most of these examples, analytical solvability relies on the absence of spatial structure, which allows for more or less conventional mean-field statistical mechanics.

In this paper we extend the class of analytically solvable models in this field. We present a model for secondary structure formation in random hetero-polymers consisting of amino-acid monomers which are allowed to interact in three qualitatively different ways: (i) via so-called steric interactions, which reflect monomer-specific geometric constraints and electrical forces determining the local energy landscape for the orientation of monomer-connecting links, (ii) via hydrogen-bonding, which acts over larger distances along the chain, and is believed to play a role in the stabilization of helix-type structures, and (iii) via polarity-induced energy gradients, which tend to promote states in which the hydrophobic amino-acids are more or less turned towards the same side of the polymeric chain, in order to enable effective shielding from water molecules via folding of the polymer as a whole. Interactions (i) and (ii) are of a short-range nature, whereas (iii) is long-range. We note that secondary structure formation has also been studied within a mean-field approach in , and that a combination of different types of monomer interactions has been proposed previously in . In the latter study, assuming statistical independence of energy levels, the random energy scheme could provide qualitative results; however, the validity of this approach has since then been questioned . In contrast, our solution does not employ random energy considerations. It is based on a combination of mean-field and random transfer-matrix techniques, which in one-dimensional models are known to reduce the evaluation of the partition function to a relatively simple numerical problem. Due to the presence of additional long-range interactions (via polarity-induced forces) our model no longer lies in the universality class of one dimensional systems, and phase transitions are therefore possible (and will indeed occur) at finite temperatures.

Our paper is organized as follows. We first define our model and the relevant macroscopic observables. Since the disordered infinite-range (polarity induced) part of our Hamiltonian, which drives the collapse to a folded state, is different in structure from the more familiar Mattis-like  terms in mean-field spin systems, we first solve our model for the case where only polarity energies are present. We then proceed to the solution of the full model, with all three interaction terms present, but now limiting ourselves (for simplicity) to the simplest choice of angular variables. Our phase diagrams exhibit second-order transitions between folded and unfolded states, whereas close to zero-temperature a hierarchy of ‘mixed’ phases appears, where new ergodic components are created and where folding depends on initial conditions. The latter phases are found to be related to entropic discontinuities. Finally, we present results from simulation experiments, which show excellent agreement with the theory.

## 2 Model Definitions

We consider one dimensional models of random hetero-polymers, where $`N`$ clock-state spin variables $`\varphi _i\{\frac{(2k+1)\pi }{q};k=0,\mathrm{},q1\}`$ describe the spatial orientations of successive monomer residues in planes vertical to the polymer’s chain axis, see figure 1. The configurational state of the system as a whole is written as $`\mathit{\varphi }=(\varphi _1,\mathrm{},\varphi _N)`$.

We define the Hamiltionian of the system to be the sum of three qualitatively different terms, $`H(\mathit{\varphi })=H_\mathrm{s}(\mathit{\varphi })+H_\mathrm{p}(\mathit{\varphi })+H_{\mathrm{Hb}}(\mathit{\varphi })`$, which are defined and interpreted as follows:

1. Polarity-induced energy (see figure 2):

$$H_\mathrm{p}(\mathit{\varphi })=\frac{J_p}{N}\underset{ij}{}\xi _i\xi _j\delta _{\varphi _i,\varphi _j}$$
(1)
This describes exchange-energies of monomer pairs generated by their polarity type, believed to be the main driving forces for compactification. Proteins live in an aqueous environment, and amino-acids of the same polarity prefer to co-align, so that folding allows the chain to arrange for hydrophobic and hydrophilic monomers to form the inner- and surface-residues of the molecule, respectively. Equation (1) describes this effect phenomenologically: $`\xi _i`$ indicates whether the monomer at site $`i`$ is hydrophobic ($`\xi _i=1`$) or hydrophilic ($`\xi _i=1`$), and we reduce the configuration energy for every pair $`(i,j)`$ of monomer residues which are both of the same type and which are also found in identical orientations relative to the backbone.
2. hydrogen-bond energy (see figure 3):

$$H_{\mathrm{Hb}}(\mathit{\varphi })=\underset{i}{}\left\{J_{Hb}^L\underset{k=0}{\overset{q1}{}}\delta _{\varphi _{i+k+1}\varphi _{i+k},\frac{2\pi }{q}}+J_{Hb}^R\underset{k=0}{\overset{q1}{}}\delta _{\varphi _{i+k+1}\varphi _{i+k},\frac{2\pi }{q}}\right\}$$
(2)
The second contribution to the energy describes the effect of hydrogen bonding: a monomer pair $`(i,j)`$ is coupled by a hydrogen bond of strength $`J_{Hb}^L`$ or $`J_{Hb}^R`$ if and only if they are spatially separated by exactly $`q`$ lattice sites and if the relative angles $`\varphi _{k+1}\varphi _k`$ of all monomers $`k=i,\mathrm{},i+q1`$ form a local helical twist of $`\pm \frac{2\pi }{q}`$ (and therefore monomer $`i`$ and monomer $`i+q`$ have the same orientation relative to the backbone), such that intermediate monomers do not block the formation of the hydrogen bond.
3. Steric energy (see figure 3):

$$H_\mathrm{s}(\mathit{\varphi })=J_s\underset{i}{}\mathrm{cos}[(\varphi _{i+1}\varphi _i)(\varphi _i\varphi _{i1})a_i]$$
(3)
This describes local short-range steric monomer-monomer interactions, favoring alignment of the relative angles $`(\varphi _{i+1}\varphi _i)`$ and $`(\varphi _i\varphi _{i1})`$ towards a specific preferred direction $`a_i`$ which depends on the type of monomer present at site $`i`$.

The various energy scales in the problem, and thus the relative importance of the three types of forces, are controlled by the non-negative coupling constants $`\{J_p,J_{Hb}^{L,R},J_s\}`$. A preference for left-or right-handed helices can be built in by modifying the balance between $`J_{Hb}^L`$ and $`J_{Hb}^R`$. The quenched disorder in the problem is given by the realisation of the (randomly drawn, but fixed) amino-acid sequence, i.e. the variables $`\{\xi _i,a_i\}`$. We denote the monomer type found at location $`i`$ in the chain by $`\lambda _i`$, so that $`\xi _i=\xi (\lambda _i)`$ and $`a_i=a(\lambda _i)`$. The disorder is characterized by the distributions

$`w[a,\xi ]`$ $`=`$ $`\underset{N\mathrm{}}{lim}{\displaystyle \frac{1}{N}}{\displaystyle \underset{i}{}}\delta _{\xi ,\xi _i}\delta [aa_i]={\displaystyle \underset{\lambda }{}}W(\lambda )\delta _{\xi ,\xi (\lambda )}\delta [aa(\lambda )]`$ (4)
$`W(\lambda )`$ $`=`$ $`\underset{N\mathrm{}}{lim}{\displaystyle \frac{1}{N}}{\displaystyle \underset{i}{}}\delta _{\lambda ,\lambda _i}`$ (5)
Note that for random hetero-polymers the distribution $`W(\lambda )`$ is simply the a priori distribution according to which the monomers were selected. The marginal distribution specifying polarity statistics is written as

$$w[\xi ]=𝑑aw[a,\xi ]=\frac{1}{2}(1+p)\delta _{\xi ,1}+\frac{1}{2}(1p)\delta _{\xi ,1},$$
(6)
with $`p[1,1]`$.

We will solve our model in thermal equilibrium via a suitable combination of mean- and random-field techniques , which will allow us to evaluate the free energy per monomer $`f`$ in the thermodynamic limit:

$$f=\underset{N\mathrm{}}{lim}\frac{1}{\beta N}\mathrm{log}\underset{\mathit{\varphi }}{}e^{\beta H(\mathit{\varphi })}$$
(7)
where $`H(\mathit{\varphi })=H_\mathrm{p}(\mathit{\varphi })+H_{\mathrm{Hb}}(\mathit{\varphi })+H_\mathrm{s}(\mathit{\varphi })`$. The parameter $`\beta `$ is an effective inverse temperature, which controls the amount of stochasticity in the underlying dynamics (with $`\beta =0`$ and $`\beta =\mathrm{}`$ corresponding to purely random and purely deterministic dynamics, respectively). The effective temperature will generally depend on various environmental factors, such as solvent conditions. We wish to emphasise that our present model takes into account the folding of the hetero-polymer only as a general mechanism with which to realise the potential for energy gain via polarity-induced forces, without specifying the detailed three-dimensional structure this reduction would give rise to. It can consequently describe only the formation of secondary structure as the result of folding, not the emerging tertiary structure; this is the price to be paid for exact analytical solvability.

Given the above definitions, it is natural to divide the monomers into two groups according to their polarity, $`\{1,\mathrm{},N\}=I_+I_{}`$ with $`I_\pm =\{i|\xi _i=\pm 1\}`$. We note that $`lim_N\mathrm{}|I_\pm |/N=\frac{1}{2}(1\pm p)`$. Within each group one can define as natural observables to measure the degree of polymer compactification (i.e. the impact of the polarity-induced forces) the distribution of monomer residue orientations, $`P_+(\varphi ;\mathit{\varphi })`$ and $`P_{}(\varphi ;\mathit{\varphi })`$:

$$P_\pm (\varphi )=\underset{N\mathrm{}}{lim}P_\pm (\varphi ;\mathit{\varphi })P_\pm (\varphi ;\mathit{\varphi })=\frac{1}{|I_\pm |}\underset{iI_\pm }{}\delta _{\varphi ,\varphi _i}$$
(8)
where $`\mathrm{}`$ denotes an average over the Boltzmann distribution $`p_{\mathrm{}}(\mathit{\varphi })\mathrm{exp}[\beta H(\mathit{\varphi })]`$. Note that, by definition, $`P_\pm (\varphi )[0,1]`$ and $`_\varphi P_\pm (\varphi )=1`$. Note also that due to the equivalence of all absolute orientations, spontaneous symmetry breaking can occur. In order to measure the degree of $`L`$(eft) or $`R`$(ight) chirality of the folded state, as induced by the steric interactions and hydrogen bonds, we introduce the two order parameters

$$\chi _\pm =\underset{N\mathrm{}}{lim}\chi _\pm (\mathit{\varphi })\chi _\pm (\mathit{\varphi })=\frac{1}{N}\underset{i}{}\underset{k=0}{\overset{m1}{}}\delta _{\varphi _{i+k+1}\varphi _{i+k},\pm \frac{2\pi }{m}}$$
(9)
Thus $`\chi _+=f/J_{Hb}^L`$ and $`\chi _{}=f/J_{Hb}^R`$.

Before solving the full model it is instructive to consider the various limiting cases one obtains by setting specific combinations of the characteristic energies $`\{J_s,J_p,J_{Hb}\}`$ in (1-3) to zero. First, in the absence of polarity interactions the model reduces to a one-dimensional random-field Potts model with site-disorder, for which the free energy is known to be analytic for finite temperatures, and there can be no phase transition. On the other hand, in the absence of steric- and hydrogen-bond interactions the model reduces to a mean-field model with site disorder. The most interesting scenario, from a physical and a technical point of view, is the one where all three forces are included. Due to the long-range interactions our model is expected to show a phase transition, whereas the short-range interactions are expected to generate frustration phenomena such as hierarchies of discontinuous transitions  and non-analytic distribution functions for local observables such as Devil’s Staircases . An appealing feature of the model is that, apart from the mean field forces, it is essentially one-dimensional and thus allows for an exact solution based on random-field techniques such as in .

## 3 Solution of the Polarity Model

In order to identify and interpret the properties of the full model, to be analysed in a subsequent section, we will now first solve our model in the absence of short-range interactions, i.e. for $`J_s=J_{Hb}^{L,R}=0`$, so that $`H(\mathit{\varphi })=H_\mathrm{p}(\mathit{\varphi })`$.

### 3.1 Calculation of the Free Energy

Upon using the simple identity $`_{ij}\delta _{\varphi _i,\varphi _j}=_\varphi _{ij}\delta _{\varphi _i,\varphi }\delta _{\varphi _j,\varphi }`$ we can express the polarity Hamiltonian (1) in terms of the order parameters (8)

$$H_\mathrm{p}(\mathit{\varphi })=J_pN\underset{\varphi }{}\left\{\frac{|I_+|}{N}P_+(\varphi ;\mathit{\varphi })\frac{|I_{}|}{N}P_{}(\varphi ;\mathit{\varphi })\right\}^2$$
(10)
Upon introducing delta functions to enforce the definitions (8), in integral representation, we obtain the following expression for the free energy per site (7):

$$f=\underset{N\mathrm{}}{lim}\frac{1}{\beta N}\mathrm{log}\underset{\varphi }{}\left[dP_\pm (\varphi )d\widehat{P}_\pm (\varphi )\right]e^{NG\left[\{P_\pm ,\widehat{P}_\pm \}\right]}$$
(11)
where

$`G\left[\{P_\pm ,\widehat{P}_\pm \}\right]`$ $`=`$ $`{\displaystyle \frac{1}{4}}\beta J_p{\displaystyle \underset{\varphi }{}}\left\{(1+p)P_+(\varphi )(1p)P_{}(\varphi )\right\}^2`$
$`i{\displaystyle \underset{\varphi }{}}\left\{\widehat{P}_+(\varphi )P_+(\varphi )+\widehat{P}_{}(\varphi )P_{}(\varphi )\right\}`$
$`{\displaystyle \frac{1}{2}}(1+p)\mathrm{log}{\displaystyle \underset{\varphi }{}}e^{2i\widehat{P}_+(\varphi )/(1+p)}{\displaystyle \frac{1}{2}}(1p)\mathrm{log}{\displaystyle \underset{\varphi }{}}e^{2i\widehat{P}_{}(\varphi )/(1p)}`$
In the thermodynamic limit $`N\mathrm{}`$, the integral in (11) can be evaluated via steepest descent. Derivation of $`G[\mathrm{}]`$ with respect to $`P_\pm (\varphi )`$ gives the equation $`i\widehat{P}_\pm (\varphi )=\frac{1}{2}(1\pm p)\beta J_p\left[(1+p)P_+(\varphi )(1p)P_{}(\varphi )\right]`$, with which we eliminate the conjugate order parameters. This results in $`f=\mathrm{extr}_{\{L\}}f[\{L\}]`$

$$f[\{L\}]=\frac{J_p}{4}\underset{\varphi }{}L^2(\varphi )\frac{1+p}{2\beta }\mathrm{log}\underset{\varphi }{}e^{\beta J_pL(\varphi )}\frac{1p}{2\beta }\mathrm{log}\underset{\varphi }{}e^{\beta J_pL(\varphi )}$$
(12)
where $`L(\varphi )=(1+p)P_+(\varphi )(1p)P_{}(\varphi )`$. Extremisation with respect to the $`L(\varphi )`$ leads to a set of $`q`$ coupled saddle-point equations from which to solve $`\{L(\varphi )\}`$, in terms of which we can then also express our original observables $`P_\pm (\varphi )`$:

$$L(\varphi )=(1+p)\frac{e^{\beta J_pL(\varphi )}}{_\varphi ^{}e^{\beta J_pL(\varphi ^{})}}(1p)\frac{e^{\beta J_pL(\varphi )}}{_\varphi ^{}e^{\beta J_pL(\varphi ^{})}}$$
(13)
$$P_\pm (\varphi )=\frac{e^{\pm \beta J_pL(\varphi )}}{_\varphi ^{}e^{\pm \beta J_pL(\varphi ^{})}}$$
(14)
Note that (13) is invariant under the transformation $`\{p,L(\varphi )\}\{p,L(\varphi )\}`$ $`\varphi `$, and that $`_\varphi L(\varphi )=2p`$.

The uniform high temperature solution, where $`L(\varphi )=L^{}=2p/q`$ for all $`\varphi `$ and therefore $`P_\pm (\varphi )=\frac{1}{q}`$ for all $`\varphi `$, always satisfies (13). Expansion of the free energy (12) around the uniform solution $`\{L^{}\}`$ allows us to determine the critical temperature $`T_c=1/\beta _c`$ where it becomes locally unstable. For perturbations $`\{\delta L\}`$ orthogonal to $`\{L^{}\}`$, i.e. for which $`_\varphi \delta L(\varphi )=0`$, we find

$$f[\{L^{}+\delta L\}]=f[\{L^{}\}]+\frac{J_p^2}{2q}(\frac{q}{2J_p}\beta )\underset{\varphi }{}\delta ^2L(\varphi )+𝒪(\delta ^3L)$$
(15)
Hence a second-order phase transition to an ordered state takes place at

$$T_c=\beta _c^1=\frac{2J_p}{q}.$$
(16)
(or at a higher temperature, as a first-order transition). This value is independent of the variable $`p=lim_N\mathrm{}\frac{1}{N}_i\xi _i`$, which measures the balance between hydrophobic and hydrophilic monomers.

Similarly we can find the system’s ground state, for any non-trivial value of $`m`$. Let us define $`L_g(\varphi )=lim_{T0}L(\varphi )`$, $`L_+=\mathrm{max}_\varphi L_g(\varphi )`$ and $`L_{}=\mathrm{min}_\varphi L_g(\varphi )`$, and let us denote the number of $`\varphi `$ for which $`L_g(\varphi )=L_+`$ as $`q_+1`$ and the number for which $`L_g(\varphi )=L_{}`$ as $`q_{}1`$ (with $`q_++q_{}q)`$. We assume $`L_{}<L_+`$; one can easily convince oneself that the alternative $`L_{}=L_+`$, i.e. the high temperature solution, will not give the ground state. Taking the $`T0`$ limit in the saddle-point equations (13) then shows that $`L_\pm =\frac{p\pm 1}{q_\pm }`$, and that $`L_g(\varphi )=0`$ for all $`\varphi `$ such that $`L_{}<L_g(\varphi )<L_+`$. Thus $`L_g(\varphi )`$ can take only one of three different values. The ground state energy per monomer, $`u=lim_{T0}f`$, can subsequently be obtained as the $`T0`$ limit of (12):

$`u`$ $`=`$ $`{\displaystyle \frac{1}{2}}J_p\underset{q_+,q_{}}{\mathrm{min}}\left\{{\displaystyle \frac{1}{2}}{\displaystyle \underset{\varphi }{}}L_g^2(\varphi )(1+p)\underset{\varphi }{\mathrm{max}}L_g(\varphi )+(1p)\underset{\varphi }{\mathrm{min}}L_g(\varphi )\right\}`$ (17)
$`=`$ $`{\displaystyle \frac{1}{4}}J_p\underset{q_+,q_{}}{\mathrm{max}}\left\{{\displaystyle \frac{(1+p)^2}{q_+}}+{\displaystyle \frac{(1p)^2}{q_{}}}\right\}={\displaystyle \frac{1}{2}}J_p(1+p^2)`$
The minimum is obtained for $`q_+=q_{}=1`$: there is one angle $`\varphi _+`$ with $`L_g(\varphi _+)=p+1`$, there is one angle $`\varphi _{}`$ with $`L_g(\varphi _{})=p1`$, and the remaining $`q2`$ orientations have $`L_g(\varphi )=0`$. The ground state, written in terms of the monomer densities $`P_\pm (\varphi )`$, is

$`P_+(\varphi _+)=1,`$ $`P_+(\varphi )=0\mathrm{for}\mathrm{all}\varphi \varphi _+`$ (18)
$`P_{}(\varphi _{})=1,`$ $`P_{}(\varphi )=0\mathrm{for}\mathrm{all}\varphi \varphi _{}`$ (19)
All hydrophobic monomers cluster at some orientation $`\varphi _+`$, and all hydrophilic monomers cluster at a different orientation $`\varphi _{}`$, which is indeed the energetically most favourable configuration for any value of $`q`$. For $`q>2`$ this introduces a trivial degeneracy of the ground state, since the choice made for $`\varphi _\pm `$ is constrained only by $`\varphi +\varphi _{}`$.

In general, non-trivial solutions of the non-linear fixed point equations (13) can only be determined numerically, due to the presence of two terms $`_\varphi e^{\pm \beta J_pL(\varphi )}`$, which act as normalisation constants for $`P_\pm (\varphi )`$ and couple the $`q`$ equations in a transcendental manner. However, for the two simplest scenarios $`q=2`$ (i.e. $`\varphi \{\frac{\pi }{2},\frac{\pi }{2}\}`$) and $`q=3`$ (i.e. $`\varphi \{\frac{2\pi }{3},0,\frac{2\pi }{3}\}`$) it turns out that these terms can be transformed away, and that an analytical solution is available. We note that, due to the specific properties of the high temperature state (where all $`L(\varphi )`$ are identical) and of the ground state (where the $`L(\varphi )`$ can take only one of three possible values), the $`q>3`$ phase diagrams can at most differ quantitatively from that of the $`q=3`$ model (provided $`q`$ remains finite). We now solve our saddle-point equations (13) for $`q\{2,3\}`$.

### 3.2 Phase Diagram for $`q=2`$

In the case where $`q=2`$ (two available orientations per monomer) we have $`\varphi \{\frac{\pi }{2},\frac{\pi }{2}\}`$, and we define $`Z=\frac{1}{2}\beta J_p[L(\frac{1}{2}\pi )L(\frac{1}{2}\pi )]`$. Since the two order parameters $`L(\varphi )`$ also obey $`\frac{1}{2}[L(\frac{1}{2}\pi )+L(\frac{1}{2}\pi )]=p`$, one simply has

$$L(\pm \frac{1}{2}\pi )=p\pm Z$$

Insertion into (13) leads to a single Curie-Weiss equation for $`Z`$:

$$Z=\mathrm{tanh}(\beta J_pZ)$$
(20)
This predicts a second-order transition at $`\beta J_p=1`$, in agreement with the critical temperature (16) for de-stabilization of the high-temperature solution found earlier. The order parameter $`Z`$ is recognised to be simply the staggered magnetisation $`N^1_i\xi _i\sigma _i`$ we would have generated if we had studied the $`q=2`$ model upon transforming $`\varphi _i=\frac{1}{2}\pi \sigma _i`$, with $`\sigma _i\{1,1\}`$ (this would have led to a Mattis-type  Hamiltonian). The order parameters $`P_+(\varphi )`$ and $`P_{}(\varphi )`$ subsequently follow in terms of the solution $`Z`$ of equation (20) as

$`P_+({\displaystyle \frac{1}{2}}\pi )=P_{}({\displaystyle \frac{1}{2}}\pi )={\displaystyle \frac{1}{1+e^{2\beta J_pZ}}}`$
$`P_+({\displaystyle \frac{1}{2}}\pi )=P_{}({\displaystyle \frac{1}{2}}\pi )={\displaystyle \frac{1}{1+e^{2\beta J_pZ}}}`$
For $`T>T_c=J_p`$ one simply recovers the uniform state $`P_+(\varphi )=P_{}(\varphi )=\frac{1}{2}`$, for all $`\varphi `$, as it should. Below $`T_c`$ the system will choose to gradually align hydrophobic and hydrophilic monomers and fold, with perfect alignment (or separation) of the two polarity types at $`T=0`$.

### 3.3 Phase Diagram for $`q=3`$

In the case where $`q=3`$ (three possible orientations per monomer) we have $`\varphi \{\frac{2}{3}\pi ,0,\frac{2}{3}\pi \}`$. The possible solutions of our saddle-point equation (13) can be classified on the basis of the number of different values taken by the three order parameters $`\{L(\frac{2}{3}\pi ),L(0),L(\frac{2}{3}\pi )\}`$, as follows:

1. All order parameters take the same value, $`L(\frac{2}{3}\pi )=L(0)=L(\frac{2}{3}\pi )=\frac{2}{3}p`$. This is the uniform high temperature state, which we have already encountered, and which according to (16) becomes locally unstable at $`T_c=\frac{2}{3}J_p`$.
2. Exactly two order parameters take the same value. In view of the invariance of equation (13) under permutations of the three allowed locations $`\{\frac{2}{3}\pi ,0,\frac{2}{3}\pi \}`$ we may without loss of generality put $`L(\pm \frac{2}{3}\pi )=L_1`$ and $`L(0)=L_2`$ (with $`L_1L_2`$).
3. All three order parameters are different: $`L(\frac{2}{3}\pi )=L_1`$, $`L(0)=L_2`$, $`L(\frac{2}{3}\pi )=L_3`$, with $`L_1L_2L_3`$.

We will show that, as the temperature is lowered, first the type (ii) solution bifurcates continuously from the type (i) solution at $`T_c^I=\frac{2}{3}J_p`$, and that the type (iii) solution, in turn, bifurcates continuously from type (ii) at a lower temperature $`T_c^{II}`$.

In order to find the type (ii) solutions, and the critical temperature for which these are created as bifurcations away from the uniform one,we introduce $`Z=L_1L_2`$. Thus

$`L(\pm {\displaystyle \frac{2}{3}}\pi )`$ $`=`$ $`L_1={\displaystyle \frac{1}{3}}(2p+Z)`$
$`L(0)`$ $`=`$ $`L_2={\displaystyle \frac{1}{3}}(2p2Z)`$
Insertion shows that such states indeed solve (13), with $`Z`$ following from

$`Z=F(Z;\beta J_p)`$ (21)
$`F(Z;K)=(1+p){\displaystyle \frac{1e^{KZ}}{2+e^{KZ}}}(1p){\displaystyle \frac{1e^{KZ}}{2+e^{KZ}}}`$
The trivial solution $`Z=0`$ of (21) brings us back to the uniform state. Bifurcations occur when $`Z=F(Z;\beta J_p)`$ and $`1=_ZF(Z;\beta J_p)`$; continuous bifurcations away from $`Z=0`$ occur when $`1=lim_{Z0}_ZF(Z;\beta J_p)=\frac{2}{3}\beta J_p`$. This gives a second-order transition from state (i) to state (ii) at the critical temperature $`T_c^I=\frac{2}{3}J_p`$, i.e. precisely at the point (16) where the uniform state was found to de-stabilise. Since $`lim_{Z\pm \mathrm{}}F(Z;K)=\pm \frac{3}{2}\frac{1}{2}p`$ and $`lim_{Z0}_Z^2F(Z;K)=\frac{2}{9}K^20`$ there is no evidence for first-order transitions.

Next, in order to analyse the type (iii) solutions and to build in the normalisation $`_\varphi L(\varphi )=2p`$, we define $`Z_1=L_1L_2`$ and $`Z_2=L_1L_3`$, such that

$`L({\displaystyle \frac{2}{3}}\pi )`$ $`=`$ $`L_1={\displaystyle \frac{1}{3}}(2p+Z_1+Z_2)`$
$`L(0)`$ $`=`$ $`L_2={\displaystyle \frac{1}{3}}(2p2Z_1+Z_2)`$
$`L({\displaystyle \frac{2}{3}}\pi )`$ $`=`$ $`L_3={\displaystyle \frac{1}{3}}(2p+Z_12Z_2)`$
This reduces our saddle-point equations (13) to two coupled equations for $`\{Z_1,Z_2\}`$, which take the following form:

$`Z_1=F(Z_1,Z_2;\beta J_p)Z_2=F(Z_2,Z_1;\beta J_p)`$ (22)
$`F(Z_1,Z_2;K)=(1+p){\displaystyle \frac{1e^{KZ_1}}{1+e^{KZ_1}+e^{KZ_2}}}(1p){\displaystyle \frac{1e^{KZ_1}}{1+e^{KZ_1}+e^{KZ_2}}}`$
For $`\{Z_1=0,Z_20\}`$ or $`\{Z_2=0,Z_10\}`$ we return to a state of type (ii), whereas the trivial solution $`Z_1=Z_2=0`$ brings us back to state (i). Bifurcations occur when $`(Z_1,Z_2)=𝑭(Z_1,Z_2;\beta J_p)`$ and $`\mathrm{det}[\mathrm{𝟏}(D𝑭)(Z_1,Z_2)]=0`$, where $`𝑭:\mathrm{}^2\mathrm{}^2`$ denotes the non-linear mapping $`(Z_1,Z_2)(F(Z_1,Z_2;\beta J_p),F(Z_2,Z_1;\beta J_p))`$ and $`D𝑭`$ its Jacobian matrix. Thus, when the system is in a type (ii) state, corresponding to e.g. $`Z_1=Z`$ and $`Z_2=0`$ with $`Z`$ given as the solution of (21), a continuous bifurcation is signaled by

$$\mathrm{det}\left|\begin{array}{cc}1(_1F)(Z,0;\beta J_p)& (_2F)(Z,0;\beta J_p)\\ (_2F)(0,Z;\beta J_p)& 1(_1F)(0,Z;\beta J_p)\end{array}\right|=0$$

Working out the partial derivatives shows that, since one of the off-diagonal terms vanishes, this is equivalent to requiring either

$$\frac{1}{\beta J_p}=\frac{1+p}{2+e^{\beta J_pZ}}+\frac{1p}{2+e^{\beta J_pZ}}$$

or

$$\frac{1}{3\beta J_p}=\frac{(1+p)e^{\beta J_pZ}}{(2+e^{\beta J_pZ})^2}+\frac{(1p)e^{\beta J_pZ}}{(2+e^{\beta J_pZ})^2}$$

The second equation signals a possible destabilisation within the class of type (ii) solutions (which can only happen when there are multiple stable type (ii) solutions, for which there is no evidence); the first equation describes the creation/annihilation of type (iii) solutions from a type (ii) one. When solved in combination with the saddle-point equation (21), this latter equation gives the desired (second-order) (ii)$``$(iii) transition line $`T_c^{II}`$. The solution can be represented conveniently in the form of a parametrisation in the $`(\beta J_p,p)`$ plane, with $`x=\beta J_pZ(\mathrm{},\mathrm{})`$:

$$\beta J(x)=\frac{1}{2}\frac{\mathrm{cosh}(x)1+x\mathrm{sinh}(x)}{\mathrm{cosh}(x)1}$$
(23)
$$p(x)=\frac{x\mathrm{cosh}(x)+2x3\mathrm{sinh}(x)}{1\mathrm{cosh}(x)x\mathrm{sinh}(x)}$$
(24)
Note that $`lim_x\mathrm{}p(x)=\pm 1`$ and that $`\beta J(x)\frac{1}{2}x`$ as $`x\mathrm{}`$. Equations (23-24) for the (ii)$``$(iii) transition, together with $`\beta J_p=3/2`$ (16) describing the (i)$``$(ii) transition, in fact represent all phase transitions in the $`q=3`$ system with polarity energies only. This conjecture is based on extensive numerical exploration of the solutions of the fixed-point equations (13).

In figure 4 we show the resultant phase diagram of the polarity model for $`q=3`$, i.e. for $`\varphi \{\frac{2\pi }{3},0,\frac{2\pi }{3}\}`$, in the $`(T/J_p,p)`$ plane. It consists of regions characterised by the number of different values taken by the three order parameters $`\{L(\varphi )\}`$, and therefore by the monomer distributions $`\{P_\pm (\varphi )\}`$ at the three possible orientations. All regions are separated by second-order transition lines, viz. (16) and (23-24). For $`T/J_p>\frac{3}{2}`$ (the high-temperature region) the only possible solution of (13), for any $`p`$, is $`L(\varphi )=\frac{2}{3}p`$ for all $`\varphi `$; here the monomers have no preferred orientation (to be interpreted as resulting in a swollen state). For $`T/J_p\frac{2}{3}`$ the equilibrium solution will depend on the value of the polarity statistics parameter $`p`$. In region (ii) the monomers exhibit some degree of orientation preference (to be interpreted as resulting in a partially folded state), whereas in region (iii) one finds a highly orientation specific solution (to be interpreted as resulting in a fully folded state). Note that, in view of the fact that also for $`q>3`$ the system will in equilibrium allow for at most three different values for the order parameters $`L(\varphi )`$, see (18,19), one must expect the $`q>3`$ phase diagrams to be qualitatively similar to the $`q=3`$ one, with only $`q`$-dependent re-scaling and weak deformations of transition lines.

In figure 5 we show the values of the three order parameters $`L(\varphi )`$, from which the monomer densities $`P_\pm (\varphi )`$ follow via (14), as a function of $`\beta J_p`$, for $`p=0.2`$ (left graph) and $`p=0.4`$ (right graph). These values are obtained by numerical solution of the saddle-point equations (13). We clearly observe the point where type (ii) solutions (two possible values for the $`L(\varphi )`$) bifurcate from the type (i) solution (all $`L(\varphi )`$ are identical), at $`\beta J_p=3/2`$ for both graphs. In contrast, the location of the second bifurcation from type (ii) to type (iii) solutions is indeed seen to depend on the parameter $`p`$, as predicted. We also observe how for $`\beta \mathrm{}`$ the system approaches the ground state (18,19), where $`L(\varphi )\{p1,0,p+1\}`$.

## 4 Solution of the Full Model for $`q=2`$

We will now turn to the full model described by the combination of all three energy contributions (1-3). Since we now have a Hamiltonian with both (site) disorder and short-range interactions, a simple mean-field approach such as that used in the previous section will no longer apply. Here our solution will be based on a suitable adaptation of the random-field techniques of . We will, for simplicity, consider only the simplest non-trivial case $`q=2`$, where $`\varphi _i=\frac{1}{2}\pi \sigma _i`$ with $`\sigma _i\{1,1\}`$. Our orientation variables can now be replaced by Ising spins, which leads to significant simplifications. For instance, the various terms in the Hamiltonian reduce to (upon dropping the irrelevant constants):

$`H_\mathrm{p}(𝝈)={\displaystyle \frac{J_p}{2N}}{\displaystyle \underset{ij}{}}\sigma _i\xi _i\xi _j\sigma _j`$ (25)
$`H_{\mathrm{Hb}}(𝝈)={\displaystyle \frac{1}{2}}J_{Hb}{\displaystyle \underset{i}{}}[1\sigma _i\sigma _{i+1}][1\sigma _i\sigma _{i1}]`$ (26)
$`H_\mathrm{s}(𝝈)=J_s{\displaystyle \underset{i}{}}\eta _i\sigma _{i+1}\sigma _{i1}`$ (27)
with $`J_{Hb}=\frac{1}{2}(J_{Hb}^L+J_{Hb}^R)`$, and with $`\eta _i=\mathrm{cos}[a_i]`$. Left- and right chirality energies have become identical, as expected for $`\varphi _i\{\frac{1}{2}\pi ,\frac{1}{2}\pi \}`$. The two ‘chirality’ order parameters (9) reduce to

$$\chi =\underset{N\mathrm{}}{lim}\chi (𝝈)\chi (𝝈)=\frac{1}{4N}\underset{i}{}[1\sigma _i\sigma _{i+1}][1\sigma _i\sigma _{i1}]$$
(28)
The joint distribution $`\stackrel{~}{w}[\eta ,\xi ]`$ of the disorder variables $`\{\eta _i,\xi _i\}`$ (which are independent for different sites) follows from (4):

$$\stackrel{~}{w}[\eta ,\xi ]=\underset{\lambda }{}W(\lambda )\delta _{\xi ,\xi (\lambda )}\delta \left[\eta \mathrm{cos}[a(\lambda )]\right]$$
(29)
### 4.1 Calculation of the Free Energy

We note that the polarity energy (25) can be written in terms of the ‘staggered magnetisation’

$$m(𝝈)=\frac{1}{N}\underset{i}{}\xi _i\sigma _i$$
(30)
in the Mattis  form $`H_\mathrm{p}(𝝈)=\frac{1}{2}J_pNm^2(𝝈)`$. We isolate the order parameter $`m`$ in the expression for the free energy per site (7), by inserting $`1=𝑑m\delta [m\frac{1}{N}_i\sigma _i\xi _i]`$. Writing the delta function in integral representation then leads to

$$f=\underset{N\mathrm{}}{lim}\frac{1}{\beta N}\mathrm{log}𝑑m𝑑\widehat{m}e^{\beta NG_N(m,\widehat{m})}$$
(31)
$$G_N(m,\widehat{m})=im\widehat{m}\frac{1}{2}J_pm^2\frac{1}{\beta N}\mathrm{log}Z_N(i\beta \widehat{m})$$
(32)
where the complicated (short-range) part of the partition sum has now been concentrated in the function $`Z_N(x)`$:

$$Z_N(x)=\underset{\sigma _1\mathrm{}\sigma _N}{}e^{\frac{1}{2}\beta J_{Hb}_i[1\sigma _i\sigma _{i+1}][1\sigma _i\sigma _{i1}]+\beta J_s_i\sigma _{i1}\eta _i\sigma _{i+1}+x_i\sigma _i\xi _i}$$
(33)
(with $`\sigma _0=\sigma _{N+1}0`$). The integral in (31) can for $`N\mathrm{}`$ be evaluated via steepest descent, and will be dominated by the saddle points of the exponent $`G_{\mathrm{}}(m,\widehat{m})`$. After elimination of $`\widehat{m}`$ via the saddle-point equation $`i\widehat{m}=J_pm`$, we can thus write the asymptotic free energy per monomer (32) as

$$f=\mathrm{extr}_m\left\{\frac{1}{2}J_pm^2\underset{N\mathrm{}}{lim}\frac{1}{\beta N}\mathrm{log}Z_N(\beta J_pm)\right\}$$
(34)
In order to calculate the partition sum (33) we will employ the random-field techniques of . We condition the function $`Z_N(x)`$ on the values $`\{\sigma _{N1},\sigma _N\}`$ of the two spins at the end of the chain:

$$Z_{\sigma \sigma ^{}}^{(N)}(x)=\underset{\sigma _1\mathrm{}\sigma _N}{}e^{\frac{1}{2}\beta J_{Hb}_i[1\sigma _i\sigma _{i+1}][1\sigma _i\sigma _{i1}]+J_s_i\sigma _{i1}\eta _i\sigma _{i+1}+x_i\sigma _i\xi _i}\delta _{\sigma _{N1},\sigma }\delta _{\sigma _N,\sigma ^{}}$$
(35)
with $`Z_N(x)=_{\sigma \sigma ^{}=\pm 1}Z_{\sigma ,\sigma ^{}}^{(N)}(x)`$. The addition of an extra monomer to the chain, i.e. $`NN+1`$, then leads to the following recurrent relation for the conditioned partition functions:

$$\left(\begin{array}{c}Z_{++}^{(N+1)}(x)\\ Z_+^{(N+1)}(x)\\ Z_+^{(N+1)}(x)\\ Z_{}^{(N+1)}(x)\end{array}\right)=𝑴_{N+1}(x)𝑻_N\left(\begin{array}{c}Z_{++}^{(N)}(x)\\ Z_+^{(N)}(x)\\ Z_+^{(N)}(x)\\ Z_{}^{(N)}(x)\end{array}\right)$$
(36)
in which the 4$`\times `$4 matrices $`𝑴_i(x)`$ and $`𝑻_i`$ are defined as

$$𝑴_i(x)=\left(\begin{array}{cccc}e^{x\xi _i}& 0& 0& 0\\ 0& e^{x\xi _i}& 0& 0\\ 0& 0& e^{x\xi _i}& 0\\ 0& 0& 0& e^{x\xi _i}\end{array}\right)$$
(37)
$$𝑻_i=\left(\begin{array}{cccc}e^{\beta J_s\eta _i}& 0& e^{\beta J_s\eta _i\beta J_{Hb}}& 0\\ e^{\beta J_s\eta _i+\beta J_{Hb}}& 0& e^{\beta J_s\eta _i+2\beta J_{Hb}}& 0\\ 0& e^{\beta J_s\eta _i+2\beta J_{Hb}}& 0& e^{\beta J_s\eta _i+\beta J_{Hb}}\\ 0& e^{\beta J_s\eta _i\beta J_{Hb}}& 0& e^{\beta J_s\eta _i}\end{array}\right)$$
(38)
As a result we can now write the short-range partition sum $`Z_N(x)`$ (33) in terms of the random matrices (37,38), where the randomness is in the $`\{\xi _i,\eta _i\}`$, as

$$Z_N(x)=\left(\begin{array}{c}1\\ 1\\ 1\\ 1\end{array}\right)\left[\underset{i=3}{\overset{N}{}}𝑴_{i+1}(x)𝑻_i\right]\left(\begin{array}{c}Z_{++}^{(2)}(x)\\ Z_+^{(2)}(x)\\ Z_+^{(2)}(x)\\ Z_{}^{(2)}(x)\end{array}\right)$$
(39)
The (random) matrix product will be evaluated in terms of the following (non-negative) stochastic quantities, which represent the different ratios of the conditioned partition sums (35):

$$k_j^{(1)}=e^{2x\xi _j}\frac{Z_{++}^{(j)}}{Z_+^{(j)}}k_j^{(2)}=e^{2x\xi _j}\frac{Z_+^{(j)}}{Z_+^{(j)}}k_j^{(3)}=e^{2x\xi _j}\frac{Z_+^{(j)}}{Z_{}^{(j)}}$$
(40)
From the recurrence relation (36) it follows that the variables $`k_j^{(\mathrm{})}`$ are, in turn, generated by iteration of the following mapping:

$$k_{j+1}^{(1)}=\frac{e^{\beta J_s\eta _j}k_j^{(1)}k_j^{(2)}+e^{\beta J_s\eta _j\beta J_{Hb}}}{e^{\beta J_s\eta _j}k_j^{(1)}k_j^{(2)}+e^{\beta J_s\eta _j+\beta J_{Hb}}}e^{\beta J_{Hb}}$$
(41)
$$k_{j+1}^{(2)}=\frac{e^{\beta J_s\eta _j}k_j^{(1)}k_j^{(2)}+e^{\beta J_s\eta _j+\beta J_{Hb}}}{e^{\beta J_s\eta _j+\beta J_{Hb}}k_j^{(2)}k_j^{(3)}+e^{\beta J_s\eta _j}}k_j^{(3)}e^{2x\xi _j}$$
(42)
$$k_{j+1}^{(3)}=\frac{e^{\beta J_s\eta _j+\beta J_{Hb}}k_j^{(2)}k_j^{(3)}+e^{\beta J_s\eta _j}}{e^{\beta J_s\eta _j\beta J_{Hb}}k_j^{(2)}k_j^{(3)}+e^{\beta J_s\eta _j}}e^{\beta J_{Hb}}$$
(43)
We now use

$$\frac{1}{\beta N}\mathrm{log}Z_N(x)=\frac{1}{\beta N}\mathrm{log}Z_{}^{(N)}(x)+𝒪(\frac{1}{N})$$
(44)
and work out the conditioned partition function $`Z_{}^{(N)}(x)`$ iteratively, via the the recurrence relation (36):

$`{\displaystyle \frac{1}{N}}\mathrm{log}Z_{}^{(N)}(x)`$ $`=`$ $`{\displaystyle \frac{1}{N}}\mathrm{log}Z_{}^{(N1)}(x){\displaystyle \frac{x\xi _N}{N}}`$ (45)
$`+{\displaystyle \frac{1}{N}}\mathrm{log}\left\{e^{\beta J_s\eta _{N1}\beta J_{Hb}}k_{N1}^{(2)}k_{N1}^{(3)}+e^{\beta J_s\eta _{N1}}\right\}`$
Further iteration of this relation gives

$$\underset{N\mathrm{}}{lim}\frac{1}{N}\mathrm{log}Z_{}^{(N)}(x)=𝑑𝒌𝑑\eta P(𝒌,\eta )\mathrm{log}\left\{e^{\beta J_s\eta \beta J_{Hb}}k^{(2)}k^{(3)}+e^{\beta J_s\eta }\right\}xp$$
(46)
with $`𝒌=(k^{(1)},k^{(2)},k^{(3)})`$, where $`p=𝑑\eta _\xi \xi \stackrel{~}{w}[\eta ,\xi ]`$ (see equation (29)), and with

$$P(𝒌,\eta )=\underset{N\mathrm{}}{lim}\frac{1}{N}\underset{i}{}\delta [\eta \eta _i]\delta [𝒌𝒌_i]$$
(47)
Provided the stochastic process (41-43) is ergodic, the distribution $`P(𝒌,\eta )`$ will be identical to the (joint) stationary distribution of the pair $`\{𝒌,\eta \}`$, i.e. we may write $`P(𝒌,\eta )=lim_N\mathrm{}\frac{1}{N}_i\delta [\eta \eta _i]\delta [𝒌𝒌_i]`$. Since $`𝒌_i`$ is always statistically independent of $`\eta _i`$ according to (41-43) ($`𝒌_i`$ depends only on those $`\eta _j`$ and $`\xi _j`$ with $`j<i`$), we have $`\delta [\eta \eta _i]\delta [𝒌𝒌_i]=\delta [\eta \eta _i]\delta [𝒌𝒌_i]`$. Hence $`P(𝒌,\eta )=P_{\mathrm{}}(𝒌|x)\stackrel{~}{w}[\eta ]`$, where $`P_{\mathrm{}}(𝒌|x)`$ is the invariant distribution of the process (41-43) (which is parametrised by $`x`$, due to the occurrence of $`x`$ in (42)) and where $`\stackrel{~}{w}[\eta ]=_\xi \stackrel{~}{w}[\eta ,\xi ]`$. We thereby find (46) being replaced by

$`\underset{N\mathrm{}}{lim}{\displaystyle \frac{1}{N}}\mathrm{log}Z_{}^{(N)}(x)`$ $`=`$ $`xp`$ (48)
$`+{\displaystyle 𝑑𝒌P_{\mathrm{}}(𝒌|x)𝑑\eta \stackrel{~}{w}[\eta ]\mathrm{log}\left\{e^{\beta J_s\eta \beta J_{Hb}}k^{(2)}k^{(3)}+e^{\beta J_s\eta }\right\}}`$
As a final consequence we can now write the free energy per monomer (34) as

$`f=\mathrm{extr}_m\{{\displaystyle \frac{1}{2}}J_pm^2+J_pmp`$ (49)
$`{\displaystyle \frac{1}{\beta }}{\displaystyle }d𝒌P_{\mathrm{}}(𝒌|\beta J_pm){\displaystyle }d\eta \stackrel{~}{w}[\eta ]\mathrm{log}[e^{\beta J_s\eta \beta J_{Hb}}k^{(2)}k^{(3)}+e^{\beta J_s\eta }]\}`$
where the invariant measure $`P_{\mathrm{}}(𝒌|x)`$ of the process (41-43) is to be solved from

$$P_{\mathrm{}}(𝒌|x)=𝑑𝒌^{}P_{\mathrm{}}(𝒌^{}|x)𝑑\eta \underset{\xi }{}\stackrel{~}{w}[\eta ,\xi ]\delta \left[𝒌(𝒌^{}|x,\eta ,\xi )\right]$$
(50)
with

$$(\begin{array}{c}_1(𝒌|x,\eta ,\xi )\\ _2(𝒌|x,\eta ,\xi )\\ _3(𝒌|x,\eta ,\xi )\end{array})=(\begin{array}{c}\frac{e^{\beta J_s\eta }k_1k_2+e^{\beta J_s\eta \beta J_{Hb}}}{e^{\beta J_s\eta }k_1k_2+e^{\beta J_s\eta +\beta J_{Hb}}}e^{\beta J_{Hb}}\\ \frac{e^{\beta J_s\eta }k_1k_2+e^{\beta J_s\eta +\beta J_{Hb}}}{e^{\beta J_s\eta +\beta J_{Hb}}k_2k_3+e^{\beta J_s\eta }}k_3e^{2x\xi }\\ \frac{e^{\beta J_s\eta +\beta J_{Hb}}k_2k_3+e^{\beta J_s\eta }}{e^{\beta J_s\eta \beta J_{Hb}}k_2k_3+e^{\beta J_s\eta }}e^{\beta J_{Hb}}\end{array})$$
(51)
In the case of the one-dimensional random-field Ising model , for which the analysis is very similar, the corresponding distribution $`P_{\mathrm{}}(𝒌)`$ is known, at least in certain parameter regions, to become highly non-trivial and acquire the form of the derivative of a Devil’s Staircase. To our knowledge, no general analytic expression has been derived to describe $`P_{\mathrm{}}(𝒌)`$ for finite temperatures. Nevertheless, for the purpose of the present paper it is only a simple numerical exercise to evaluate $`P_{\mathrm{}}(𝒌|x)`$ directly by iteration of (51), for values of $`\{\eta ,\xi \}`$ drawn randomly according to $`\stackrel{~}{w}[\eta ,\xi ]`$.

### 4.2 Simple Limiting Cases

Before converting our general results into phase diagrams we will first carry out benchmark tests of our expressions, by inspecting simple limits.

* Firstly, in the absence of short-range interactions, i.e. for $`J_{Hb}=J_s=0`$, the expression for the asymptotic free energy per site (49) should reduce to the $`q=2`$ version of (12), which ought to be simply the free energy of the infinite-range Mattis model . Indeed, we find that for $`J_s=J_{Hb}=0`$ the mapping (51) reduces to

$$(\begin{array}{c}_1(𝒌|x,\eta ,\xi )\\ _2(𝒌|x,\eta ,\xi )\\ _3(𝒌|x,\eta ,\xi )\end{array})=(\begin{array}{c}1\\ \frac{k_1k_2+1}{k_2k_3+1}k_3e^{2x\xi }\\ 1\end{array})$$
(52)
Hence $`P_{\mathrm{}}(𝒌|x)=\delta [k_11]\delta [k_31]P_{\mathrm{}}(k_2|x)`$, with

$$P_{\mathrm{}}(k_2|x)=\frac{1}{2}(1+p)\delta [k_2e^{2x}]+\frac{1}{2}(1p)\delta [k_2e^{2x}]$$
(53)
Substitution into (49), for $`J_{Hb}=J_s=0`$, gives

$`f=\mathrm{extr}_m\left\{{\displaystyle \frac{1}{2}}J_pm^2{\displaystyle \frac{1}{\beta }}\mathrm{log}2\mathrm{cosh}(\beta J_pm)\right\}`$ (54)
which is indeed the well-known asymptotic free energy per site of an infinite-range Mattis magnet .
* Secondly, for $`J_{Hb}=0`$ and $`\eta _i=\xi _i=1`$ for all $`i`$ (i.e. $`\stackrel{~}{w}[\eta ,\xi ]=\delta [\eta 1]\delta [\xi 1]`$ and $`p=1`$) the macroscopic laws of our model should reduce to those of the , which describes pattern recall in recurrent neural networks with competition between short-range and long-range information processing, for the simplest ‘one-pattern’ scenario. For $`J_{Hb}=0`$ and $`\stackrel{~}{w}[\eta ,\xi ]=\delta [\eta 1]\delta [\xi 1]`$ the mapping (51) becomes fully deterministic, and takes the form

$$(\begin{array}{c}_1(𝒌|x)\\ _2(𝒌|x)\\ _3(𝒌|x)\end{array})=(\begin{array}{c}\frac{e^{\beta J_s}k_1k_2+e^{\beta J_s}}{e^{\beta J_s}k_1k_2+e^{\beta J_s}}\\ \frac{e^{\beta J_s}k_1k_2+e^{\beta J_s}}{e^{\beta J_s}k_2k_3+e^{\beta J_s}}k_3e^{2x}\\ \frac{e^{\beta J_s}k_2k_3+e^{\beta J_s}}{e^{\beta J_s}k_2k_3+e^{\beta J_s}}\end{array})$$
(55)
and $`P_{\mathrm{}}(𝒌|x)=\delta [𝒌𝒌^{}(x)]`$, where $`𝒌^{}(x)`$ denotes the fixed-point of the mapping (55) with non-negative components, which (in line with our previous assumption of ergodicity of the original process (41-43)) we assume to be unique. We observe that (55) preserves $`k_1=k_3`$, and the remaining components of the fixed-point $`𝒌^{}(x)=(k_1^{},k_2^{},k_1^{})`$ must obey

$$(\begin{array}{c}k_1^{}\\ k_2^{}\end{array})=(\begin{array}{c}\frac{e^{\beta J_s}k_1^{}k_2^{}+e^{\beta J_s}}{e^{\beta J_s}k_1^{}k_2^{}+e^{\beta J_s}}\\ \frac{e^{\beta J_s}k_1^{}k_2^{}+e^{\beta J_s}}{e^{\beta J_s}k_2^{}k_1^{}+e^{\beta J_s}}k_1^{}e^{2x}\end{array})$$
(56)
This (in turn) gives $`(k_1^{},k_2^{})=(k^{},e^{2x})`$, where (upon substituting $`x=\beta J_pm`$) $`k^{}`$ is the non-negative solution of

$$k^{}=\frac{e^{\beta (J_s+J_pm)}k^{}+e^{\beta (J_s+J_pm)}}{e^{\beta (J_sJ_pm)}k^{}+e^{\beta (J_sJ_pm)}}$$
(57)
Insertion into (49) gives us

$$f=\mathrm{extr}_m\left\{\frac{1}{2}J_pm^2\frac{1}{\beta }\mathrm{log}\left[e^{\beta (J_sJ_pm)}k^{}+e^{\beta (J_sJ_pm)}\right]\right\}$$
(58)
It follows from (57) that the quantity $`\lambda =e^{\beta (J_sJ_pm)}k^{}+e^{\beta (J_sJ_pm)}`$ occurring in (58) obeys $`(\lambda e^{\beta (J_2+J_pm)})(\lambda e^{\beta (J_sJ_pm)})=e^{2\beta J_s}`$, which we recognise as the eigenvalue equation of the transfer matrix

$$𝑻=\left(\begin{array}{cc}e^{\beta (J_s+J_pm)}& e^{\beta J_s}\\ e^{\beta J_s}& e^{\beta (J_sJ_pm)}\end{array}\right)$$
(59)
This shows that the free energy (58) is indeed identical to that of .

### 4.3 Phase Diagrams and Comparison with Numerical Experiments

In order to obtain phase diagrams we finally have to calculate the local extrema of a free energy surface $`f[m]`$, the argument of the extremisation in (49), which still depends on the choice made for the statistics of the monomer properties $`\{\xi ,\eta \}`$. Here we apply our theory to the simple example $`\stackrel{~}{w}[\eta ,\xi ]=\frac{1}{4}[\delta (\eta +1)+\delta (\eta 1)][\delta (\xi +1)+\delta (\xi 1)]`$, hence also $`p=0`$. In this case the free energy surface $`f[m]`$ simplifies to

$`f[m]`$ $`=`$ $`{\displaystyle \frac{1}{2}}J_pm^2{\displaystyle \frac{1}{2\beta }}{\displaystyle 𝑑𝒌P_{\mathrm{}}(𝒌|\beta J_pm)}`$ (60)
$`\times \mathrm{log}\left[(e^{\beta (J_s+J_{Hb})}k_2k_3+e^{\beta J_s})(e^{\beta (J_sJ_{Hb})}k_2k_3+e^{\beta J_s})\right]`$
where the invariant measure $`P_{\mathrm{}}(𝒌|x)`$ of the process (41-43) is to be solved from

$$P_{\mathrm{}}(𝒌|x)=\frac{1}{4}𝑑𝒌^{}P_{\mathrm{}}(𝒌^{}|x)\underset{\eta =\pm 1}{}\underset{\xi =\pm 1}{}\delta \left[𝒌(𝒌^{}|x,\eta ,\xi )\right]$$
(61)
with the mapping defined in (51). We determine the solution of (61) via numerical iteration. Note that, due to $`\stackrel{~}{w}[\xi ]=\stackrel{~}{w}[\xi ]`$, we have $`P_{\mathrm{}}(𝒌|x)=P_{\mathrm{}}(𝒌|x)`$. Hence $`f[m]=f[m]`$, and $`m=0`$ always corresponds to a saddle-point of $`f[m]`$. Note also that for $`J_{Hb}=0`$ (no hydrogen bonding) considerable further simplification of (60,61) will be possible, due to the resulting conservation of the symmetry $`k_1=k_3`$ by the map (51).

Examples of the results of our analysis of the surface (60) are shown in figure 6, as phase diagram cross-sections in the $`(T,J_p)`$ plane, for $`\{J_s=4,J_{Hb}=0\}`$ (left picture) and $`\{J_s=4,J_{Hb}=\frac{1}{2}\}`$ (right picture). They involve

a high-temperature phase ‘P’, where $`m=0`$ is the only local minimum of $`f[m]`$ and no folding will occur,

a phase ‘F’ where two equivalent $`m0`$ solutions minimise $`f[m]`$ (one positive, one negative, reflecting the symmetry of the present model under overall reflection $`\varphi _i\varphi _i+\pi `$), the ‘folded state’,

phases ‘M’ where four $`m0`$ solutions minimise $`f[m]`$ locally (two positive, two negative).

In the M phases, the degree of folding observed will strongly depend on initial conditions (in spite of the fact that the lowest value for $`f[m]`$, and hence the thermodynamic state, corresponds only to the maximally folded state, where $`|m|`$ is largest). See also figure 9 below. The P$``$F transition is an ordinary second-order transition, whereas the F$``$M transitions are first-order (dynamical) transitions. In the presence of hydrogen bonding, the M phases are found to be increasingly suppressed (see right picture).

In order to illuminate the physical mechanism which produces the ‘mixed’ phases, we plot in figure 7 the entropy per monomer $`s=f/T`$ close to $`T=0`$, for each of the local minima of $`f[m]`$. It is seen to become non-zero, and to develop a hierarchy of sharp peaks as a function of $`J_p`$ (c.f. ). These peaks correspond to special parameter values for which frustration effects become dominant, and for which many energetically equivalent states are possible. The largest value of the ground state entropy is obtained at the first of these peaks, for $`J_p11.2`$; this corresponds to the location in the phase diagram where the first of the ‘mixed’ phases appears, see figure 6.

The qualitative features of diagrams such as those shown in figure 6 can now be understood as follows. For large values of $`\{J_p,T\}`$ the short-range forces (steric forces and hydrogen-bonds) become irrelevant, and the diagram approaches that of a Mattis model (as it should), with a second-order transition along the line $`T=J_p`$. For low temperatures the simple Mattis state is disrupted by the steric interactions, which try to enforce monomer-specific short-range order along the chain; as a result the value needed for $`J_p`$ to create $`m0`$ states is increased (explaining the re-entrance observed in figures 6. The complex phenomenology (reminiscent of random field models) of multiple locally stable configurations, induced by the steric interactions, is subsequently found to be damped by the hydrogen bonds, which act to reduce the complexity of the ground state.

Next, in figure 8 we plot the equilibrium values of the ‘chirality’ (28) and ‘polarity’ (30) order parameters as functions of the hydrogen bond strength $`J_{Hb}`$, for three different values of $`J_p`$ (in a region of the phase diagram where there are no mixture phases, i.e. where apart from overall reflection, the stationary state is unique). Note that $`\chi `$ is simply calculated as $`\chi =\frac{1}{2}f/J_{Hb}`$ (which is done numerically). The two order parameters $`\chi `$ and $`m`$ are seen to show an opposite dependence on $`J_{Hb}`$ (monotonically increasing vs. decreasing), as they should, since, $`\chi `$ measures the degree of helical structure along the chain, whereas $`m`$ measures the probability to find monomers with identical polarity at the same side of the chain. Due to the competing roles played by two coupling parameters $`\{J_p,J_{Hb}\}`$, we see that ‘helices’ are favoured for large $`J_{Hb}`$ or small $`J_p`$ whereas ‘folding’ in the sense of efficient polarity separation, on the other hand, is favoured for small $`J_{Hb}`$ or large $`J_p`$. Note that the observed incompatibility of helical structure with polarity separation is just a reflection of the simple form we chose in this section for the disorder distribution $`\stackrel{~}{w}[\eta ,\xi ]`$ (with statistically independent $`\eta `$ and $`\xi `$); the situation would obviously have been different for distributions describing correlated disorder variables. In the same figures we also show the results of numerical simulations, for comparison (the markers in the two graphs). For small $`J_{Hb}`$ our experiments are seen to be in excellent agreement with the theory (finite size effects are of the order of $`𝒪(N^{1/2})0.03`$) whereas for large $`J_{Hb}`$ short-range couplings become increasingly dominant, leading to domain formation and very slow equilibration times, which make it difficult in practice to probe the equilibrium regime. In our experiments we have measured the value of the order parameters after $`120,000`$ iterations per spin, which for large $`J_{Hb}`$ is no longer sufficient. Note that the theory also predicts the existence of repeated small discontinuous in both order parameters; these originate from frustration-related short-range phenomena, as described in e.g. , which induce discretisation of observable supports  and non-analytic integrated distribution functions (e.g. the Devil’s Staircase ).

To verify our results further we have also performed simulation experiments in the ‘mixed’ phase regions, where our theory predicts that the extent of polarity-driven folding (i.e. the equilibrium value of $`m`$) will depend on initial conditions. In figure 9 we show the value of the ‘polarity’ order parameter $`m`$, as measured in numerical simulations of an $`N=1000`$ chain after $`20,000`$ iterations per monomer, as a function of its initial value $`m(t=0)`$, for two different parameter settings (one, to the left, in an M region of the phase diagram; one, to the right, in an F region of the phase diagram). In the insets of these graphs we also plot the corresponding free energy per monomer $`f[m]`$ as predicted by our theory, which shows either two $`m>0`$ locally stable states (left picture) or one $`m>0`$ locally stable state (right picture), respectively. In both cases the numerical experiments are found to verify the existence and the quantitative properties of the expected ergodicity breaking in the M phase. We clearly observe that, in phase M, the choice of initial conditions, in particular whether or not $`m(t=0)`$ is to the left of the free energy barrier in $`f[m]`$, determine the equilibrium value of $`m`$. We also see that in the ‘mixed’ phase (left picture) the ergodic component with the smallest value of $`m`$ is poorly equilibrated due to domain formation. This has also been observed for a similar type of statistical mechanical model in : in those parameter regions where a multiple number of states can be locally stable, different ergodic components are found to have different equilibration time-scales.

## 5 Discussion

In this paper we have presented an exactly solvable model for secondary structure formation in random hetero-polymers, consisting of amino-acid monomers which are allowed to interact in three qualitatively different ways: via (short-range) steric interactions, via (short-range) hydrogen-bonding, and via (long-range) polarity-induced forces. Our strategy was to exploit the one-dimensional nature of the monomer chain, and to separate questions relating to secondary structure formation from those relating to tertiary structure formation by taking into account the effects of the latter only via an effective energy term which measures the potential for overall energy reduction by folding (rather than trying to find the actual state realising this potential). This allows us to move away from real-space calculations towards a calculation in $`1+\mathrm{}`$ dimensions, where the statistical mechanical variables represent the orientations of the monomer residues relative to the chain axis. Solution can now be based on a combination of mean-field and random transfer-matrix techniques, which in one-dimensional models are known to reduce the evaluation of the partition function to a relatively simple numerical problem. Due to the presence of long-range interactions (via polarity-induced forces), phase transitions are still possible (and do indeed occur) at finite temperatures.

Our order parameters measure the degree of polarity-induced collapse of the chain, as well as the degree of helicity along the chain. The phase diagrams exhibit second-order transitions between ‘folded’ and ‘unfolded’ states, and, for low temperature and sufficiently strong steric interactions, a series of ‘mixed’ phases (separated from the previous ones by discontinuous transitions) where, in addition to the maximally folded states, specific partially folded states can also be locally stable. The latter phases are created at parameter values for which frustration is maximal, and where the entropy becomes particularly large. Although in the present paper we have mostly restricted ourselves (for simplicity) to chains with just a small number of possible orientations per monomer, it is not fundamentally more difficult to solve the model for larger degrees of orientational freedom (although certain adaptations are needed before the continuum limit can be taken, such as a re-scaling of the effective long range coupling $`J_p`$ and/or of the number of relative monomer orientations where polarity interactions occur). We have only evaluated our theory for the simplest choice of disorder statistics (the statistical properties of the monomers, and their physical properties such as polarity and steric constraints). Here the emerging picture is already quite satisfactory, in that explicit analytical results can be obtained, and that the predicted physical behaviour of the monomer chain (confirmed qualitatively and quantitatively by numerical simulations) makes perfect sense in the context of proteins: the polarity forces drive the transition to a collapsed state, the steric forces introduce monomer specificity, and the hydrogen bonds stabilise the conformation by damping the frustration-induced multiplicity of states.

There is still much scope for increasing the biological realism and relevance of our model without affecting its analytical solvability, at different levels. Firstly, without changing the model or its techniques for solution, one can easily consider more realistic choices for the monomer statistics, such us non-binary polarity variables, or for the orientational freedom of the monomers (for instance, the hydrogen-bond term may be modified to favour helix-type formations at the biologically observed ratio of 3.6 monomers per turn). Secondly, at a next level of sophistication one could construct a more realistic form for the polarity induced energy contribution (breaking the present hydrophobic-hydrophilic symmetry, and based upon biological data), or more realistic representations of the degrees of freedom of the individual peptide units and residues (i.e. three angles per monomer, rather than one), or the action of ‘chaperones’ (via external fields). Solution of such models would not be essentially more difficult than that of the examples worked out here; the main problem would rather be to extract the canonical definitions of the ingredients to be incorporated into the model from the available biological data. In contrast, qualitatively different and more difficult types of modification and extension would be to consider non-random hetero-polymers, where the monomer properties and statistics are chosen such as to mimic real proteins, or to try to analyse the interplay between secondary and tertiary structure formation. Here new techniques for solution will have to come in.

The main problem in the statistical mechanical study of folding proteins appears to be the construction of models where an acceptable and productive balance can be found between analytical solvability and biological realism. We believe that our present model might point to a new direction where this might be achieved.

## References