On phase diagram and the pseudogap state in a linear chiral homopolymer model A. Sinelnikova,1, ∗ A.J. Niemi,2, 3, 4, † and M. Ulybyshev5, 1, 6, ‡ 1 arXiv:1501.05267v1 [cond-mat.soft] 21 Jan 2015 Institute for Theoretical Problems of Microphysics, Moscow State University, Moscow, 119899 Russia 2 Department of Physics and Astronomy, Uppsala University, P.O. Box 803, S-75108, Uppsala, Sweden 3 Laboratoire de Mathematiques et Physique Theorique CNRS UMR 6083, F´ed´eration Denis Poisson, Universit´e de Tours, Parc de Grandmont, F37200, Tours, France 4 Department of Physics, Beijing Institute of Technology, Haidian District, Beijing 100081, P. R. China 5 Institute of Theoretical Physics, University of Regensburg, D-93053 Germany, Regensburg, Universitatsstrasse 31 6 ITEP, B. Cheremushkinskaya str. 25, Moscow, 117218 Russia The phase structure of a homopolymer chain is investigated in terms of a universal theoretical model, designed to describe the infrared limit of slow spatial variations. The effects of chirality are studied and compared with the influence of a short-range attractive interaction between monomers, at various ambient temperature values. In the high temperature limit the homopolymer chain is in the self-avoiding random walk phase. But at low temperatures, two different phases are possible: When short-range attractive interactions dominate over chirality, the chain collapses into a spacefilling conformation. But when the attractive interactions become weaker, there is a low temperature unfolding transition and the chain becomes like a straight rod. Between the high temperature and low temperature limits, several intermediate states are observed. For sufficiently high values of short-range attraction , the conventional θ-regime is observed between the self-avoiding random walk phase and the space filling collapsed phase. But when chirality increases, there is a trasition from the θ-regime to a pseudogap state. Moreover, a regime akin the θ-regime is identified between the pseudogap state, and the low temperature phase where the chain is like a straight rod. Applications to polymers and proteins, in particular collagen, are suggested. I. INTRODUCTION A linear homopolymer is made of a single type of a repeat unit. An important example is polyacetylene, an organic conductive polymer which is the paradigm material for fractional fermion number [1, 2]. Additional examples, among many others, are poly-L-lysine and polyL-glutamic. The former is a food preservative with potential for wider, even pharmaceutical antimicrobial effects [3] while the latter is used for drug delivery against cancer [4]. From a theoretical point of view, the concept of a homopolymer chain can provide a useful coarse grained approximation, even in the case of a heteropolymer that exhibits only approximatively repeating patterns: For a sufficiently long chain the distinct monomers are simply combined into appropriate subunits, to dispose of the inhomogeneities in the monomer species. For example collagen, which is the most abundant protein in mammals, displays a repeated glycine-proline-X pattern, where X is any amino acid other than glycine and proline. DNA, RNA and the Cα backbone of a protein chain are additional examples where a homopolymer approximation is occasionally profitably introduced [5–7]. Here the phase structure of linear homopolymers is in- ∗ Electronic address: ann.sinelnikova@gmail.com address: Antti.Niemi@physics.uu.se; URL: http:// www.folding-protein.org ‡ Electronic address: ulybyshev@goa.bog.msu.ru † Electronic vestigated, in terms of a universal energy function [8– 13]. Of particular interest is the effect of chirality on the phase diagram. A chiral hopolymer is one where parity is broken; the mirror image of a stable chiral homopolymer conformation is in general not stable. Chiral polymers often have a tendency to form helical structures. For example in the case of proteins right-handed helical structures are more common than left-handed ones. Even though there are also chiral proteins that can form different right-handed and left-handed structures. An important example is collagen where the right-handed polyproline I conformation is more compact than the left-handed polyproline II. A general argument which is due to Kadanoff and Wilson [14–18] states, that in the thermodynamical limit the phase transition properties of a material system are often universal i.e. independent of the atomic level details. From this perspective, the construction of the phase diagram of a linear and chiral, structureless homopolymer should be relevant for the understanding of the phase diagram of more elaborate linear chiral homopolymers, maybe even that of certain heteropolymers [5–7]. All linear polymers are presumed to have a very similar phase structure, quite independently of their chemical composition [5–7]. Even though the phase where a particular polymer resides depends on many factors such as concentration, the quality of solvent, ambient temperature and pressure. Three different, universal phases are commonly identified. These phases are categorised by the way how the polymer structure fills the space [5–7]: Under poor solvent conditions or at low temperatures, when the attractive interactions between the monomers 2 dominate, a single polymer chain is presumed to collapse into a configuration which is space filling. On the other hand, in a good solvent or at high temperatures, when the repulsive interactions dominate and cause the chain to effectively swell, its geometric structure resembles that of a self-avoiding random walk (SARW). Between the two, there is a θ-regime where the attractive and repulsive interactions cancel each other. In the θ-regime the polymer chain is presumed to have the characteristics of an ordinary random walk (RW). Finally, some polymers such as collagen for example, are more like straight, riding rods. Each of the four phases – rigid rod, SARW, RW and the space filling one – is characterised by the inverse of their Hausdorff dimension, called the scaling exponent ν [5–7]. To define this quantity, one introduces the radius of gyration Rgyr . With ri the coordinates of the individual monomers, the radius of gyration is [5–7], s 1 X Rgyr = (ri − rj )2 (1) 2N 2 i,j In the limit where the number of monomers N becomes very large, the radius of gyration has the asymptotic expansion [19–22] N large R02 N 2ν (1 + R1 N −δ1 + ...) ∼ R02 N 2ν + . . . (2) Here the length scale R0 is called the Kuhn length. It is the effective distance between the monomers in the large-N limit. The Kuhn length is not a universal quantity, its value can in principle be computed from the atomic level details of the polymer and environment including pressure, temperature and chemical microstructure of the solvent. The dimensionless scaling exponent ν i.e. the inverse Hausdorff dimension that governs the large-N asymptotic form of equation (2), is presumed to be a universal quantity. Its numerical value is independent of the local atomic level structure of the polymer [5, 19? –21]. The δ1 etc. are critical exponents and the R1 etc. are the corresponding amplitudes, and together they constitute the finite-size corrections. The δ1 etc. are universal quantities [21], but the R1 etc. are not universal [21]. The following mean field values are conventionally assigned to ν [5, 7]: 2 Rgyr −→ 1/3 1/2 ν = 3/5 1 collapsed SARW RW rod (3) Under poor solvent conditions or at low temperatures, the polymer collapses into a space-filling conformation [23–25] with the mean field exponent ν = 1/3. Folded proteins are commonly found in this phase. For an ordinary random-walk (RW) the mean field value is ν = 1/2. This corresponds to the θ regime, that separates the collapsed phase from the high-temperature self-avoiding random walk phase for which the Flory value ν = 3/5 is found. Finally, when ν = 1, the polymer loses its inherently fractal structure and behaves like a straight rod [25]. The transition between the collapsed phase and the SARW phase can involve the RW phase as a tri-critical θ-point. More generally, the RW phase can appear as a transitional θ-regime, separating the high temperature SARW and low temperature collapsed phases. The transitions between the rigid rod phase and the other three, is less studied. However, there is a physically and biologically very important scenario where such a transition could have a rˆole; that of cold denaturation of a protein chain. The presence of all four phases (3) opens the possibility of a 4-critical point, under proper conditions [26]. The θ-point value ν = 1/2 is exact for a polymer with no long range interactions [5]. For a space-filling structure the value ν = 1/3 is also exact, and similarly ν = 1 is exact for a straight, linear rod-like structure. But in the case of SARW the mean field value is corrected by fluctuations. Using the argument that the SARW is in the same universality class with the O(n) symmetric scalar field theory with a quartic self-interaction, the following modified value was obtained [19, 20] ν = 0.5880 ± 0.0015 A subsequent numerical Monte Carlo evaluation, computed directly by using the self-avoiding random-walk model on a square lattice, gave a very similar value [21] ν = 0.5877 ± 0.0006 In the present article, the phase structure of a homopolymer is studied using numerical simulations, in the context of a universal off-lattice energy function. Different Monte Carlo algorithms are introduced and analysed, to identify an optimal one for describing a homopolymer. The phase digram that corresponds to the geometric characteristics (3) is then constructed, by inspecting the radius of gyration in terms of ambient temperature and the various parameters that appear in the Hamiltonian, in particular the one that characterises the chirality and the one that characterises the attractive (hydrophobic) forces. II. THE MODEL A. Geometry Consider a homopolymer chain with i = 1, ..., N monomers and let ri be the three dimensional space coordinates of the monomers. The unit tangent vectors along the lines that connect two consecutive monomers are ti = ri+1 − ri |ri+1 − ri | (4) 3 τi → τi + ∆i−1 − ∆i The unit binormal vectors are defined by bi = ti−1 × ti |ti−1 × ti | (5) The unit normal vectors are defined by ni = bi × ti (6) The three vectors (ni , bi , ti ) determine an orthonormal frame at the monomer position ri . The discrete bond angles are κi ≡ κi+1,i = arccos (ti+1 · ti ) (7) and the discrete torsion angles are (14) The range of τi is [−π, π) mod(2π). The equations (13) and (14) may be used to extend the range of the bond angle from [0, π) to κi into [−π, π) mod(2π). The extension is compensated for by the following discrete Z2 symmetry κk → − κk τi → τi − π for all k ≥ i (15) that leaves the chain intact. In the numerical simulations presented here, all the distances between nearest neighbour monomers are fixed to the uniform constant value τi ≡ τi+1,i = sgn[(bi−1 ×bi )·ti ]×arccos (bi+1 · bi ) (8) |ri+1 − ri | = δ = 3.8 (˚ A) Conversely, when the angles (κi , τi ) are known the discrete Frenet equation [27] ni+1 cos κ cos τ cos κ sin τ − sin κ ni bi+1 = − sin τ bi cos τ 0 ti+1 sin κ cos τ sin κ sin τ cos κ i+1,i ti (9) determines the frames iteratively, by computing the frame at the position of the (i + 1)th monomer from the frame at the position of the ith monomer. Once all the frames have been constructed, the entire chain is obtained as follows, This equals the average distance between two consecutive Cα atoms along a protein backbone, measured in ˚ Angstr¨om’s. A polymer is subject to steric constraints, due to overlapping electron clouds and various short range Born repulsions. Accordingly the following forbidden volume constraint is introduced, rk = k−1 X |ri+1 − ri | · ti |ri − rk | ≥ δ ≡ 3.8 (˚ A) for |i − k| ≥ 2 (16) (17) This is in line with the minimum distance observed between any two Cα atoms, in folded protein structures. The numerical values (16) and (17) can both be independently modified, with no effect to conclusions. (10) i=0 With no loss of generality one can set r0 = 0, and orient t0 to point into the direction of the positive z-axis. A framing is necessary for the construction of the chain from the bond and torsion angles. But the equation (10) does not involve the vectors ni and bi . Thus any linear combination of these two vectors could be chosen to define a framing, to construct the chain from the angles. For this, consider a local SO(2) transformation that rotates the frame (ni , bi ) by an angle ∆i leaving ti intact, n n cos ∆ sin ∆ 0 n i i 3 b → e∆i T b = − sin ∆i cos ∆i 0b (11) t i t i 0 0 1 t i where (T a )bc = abc are the SO(3) generators [T a , T b ] = abc T c When the vectors ni and bi are combined into a complex vector, (11) can be recast as follows [27], ni + ibi → e i∆i (ni + ibi ) ≡ e1i + ie2i (12) The frame rotation (11) introduces transformations in the equation (9), that can be summarised as follows: 3 κi T 2 → e∆i T (κi T 2 ) e−∆i T 3 (13) B. Hamiltonian The bond and torsion angles constitute a complete set of geometric variables in the case of protein Cα backbones [28]. Furthermore, according to (2) (3) the structural phase diagram of a homopolymer is determined by the three dimensional chain geometry. Thus the bond and torsion angles are a complete set of order parameters, in the sense of Kadanoff and Wilson: The geometrically defined, structural phase diagram of a homopolymer can be fully determined by a thermodynamical free energy which is constructed from these order parameters only. Consider a homopolymer chain in thermal equilibrium and let F be the thermodynamical Helmholtz free energy. Thus, the minimum energy configuration of F describes the chain, under thermodynamical equilibrium conditions. The free energy is the sum of the internal energy U and the entropy S, at temperature T F = U − TS (18) It is a function of all the inter-atomic distances F = F (rαβ ) ; rαβ = |rα − rβ | (19) where the indices α, β, ... extend over all the atoms in the homopolymer system, including those of the solvent 4 environment. Consider the infrared, long distance limit where the characteristic length scales of spatial deformations along the homopolymer chain around its thermal equilibrium configuration are large in comparison to the distance (16) between neighboring monomers. This is synonymous to an assumption that there are no abrupt wrenches and buckles along the chain, that there are only gradual long wavelength bends and twists (which is the limit of adiabatic deformations). The completeness of the bond and torsion angles implies that, in order to determine the thermodynamical phase state of the homopolymer chain, it is sufficient to consider the response of all the distances between all the atoms to the variations in the bond and torsion angles only, rαβ = rαβ (κ, τ ) Here, and in the sequel, (κ, τ ) denotes collectively all the variables κi and τi . Suppose that at a local extremum of the free energy, the bond and torsion angles along the homopolymer chain have the values The expansion (21) may be re-arranged in terms of of the differences in the angles, as follows: X F (ρ) = { Vk (ρk ; ρ0k ) + k + Zk (ρk ; ρ0k )(ρk ρk+1 + ρk ρk−1 ) + . . . } (22) Here ρ0k denotes a combination of the various parameters (κi0 , τi0 ) along the chain. But Vk (ρk ; ρ0k ), Zk (ρk ; ρ0k ) and so forth depend on the variable ρk only on the site k; these functions are ultralocal. The terms that are not shown explicitly, consist of higher order differences ρk ρk+i with i ≥ 2, and higher powers of the differences. The local terms Vk (ρk ) constitute the effective potential X Vef f = Vk (ρk ) (23) k Start by Taylor expanding the infrared limit Helmholtz free energy (18) around the extremum, The structure of the effective potential is commonly used to conclude whether a spontaneous symmetry breaking takes place. The transition from (21) to (22) involves, a priori, an infinite re-arrangement of the terms in the Taylor expansion (21). In particular, the expansion (22) has been designed so that in the continuum limit where distance between neighboring monomers vanishes i.e. δ → 0 in (16), it becomes, at least naively, an expansion of the free energy in powers of momentum about the point where momentum vanishes: For a single scalar variable ρk with continuum limit F [rαβ = rαβ (κi , τi )] ≡ F (κ, τ ) = F (κ0 , τ0 ) ρk → φ(x) X ∂F ∂F ∆κk + ∆τk + ∂κk |0 ∂τk |0 the corresponding continuum limit of (22) is the derivative expansion [29] Z 1 F (φ) = V (φ) + (∂µ φ)2 Z(φ) + . . . 2 (κi , τi ) = (κi0 , τi0 ) Consider a conformation where the (κi , τi ) deviate from these extremum values. The deviations are ∆κi = κi − κi0 ∆τi = τi − τi0 (20) k + X 1 ∂2F ∂2F ∆κk ∆κl + ∆κk ∆τl + 2 ∂κk ∂κl |0 ∂κk τl |0 The free energy must remain invariant under the local frame rotations (13), (14); the physical properties of a chain can not depend on the framing. Accordingly, it has been concluded [8–13] that to the leading non-trivial order the free energy must have the form k,l 1 ∂2F + ∆τk ∆τl 2 ∂τk ∂τl |0 + O(∆3 ) H=− The first term in the expansion evaluates the free energy at the extremum. Since (κi0 , τi0 ) correspond to the extremum, the second term vanishes. Denote in the sequel (κi , τi ) collectively, as the variable ρi . Then, = F (ρ0 ) + 1 X ∂2F ∆ρk ∆ρl + O(∆3 ) 2 ∂ρk ∂ρl |0 k,l 2κi+1 κi i=1 + F (κ, τ ) ≡ F (ρ) N −1 X i=N X 2 2κi + q(κ2i − m2 )2 + cdκ2i τi2 i=1 (21) + i=N X i=1 2 cτi − aτi − bκ2i τi (24) 5 In (24) q, m, a, b, c, d depend on the atomic level physical properties and the chemical microstructure of the homopolymer chain and its environment. In principle, these parameters can be computed from this knowledge. (Note the combination cd in the last term of second sum; this choice is made for later convenience.) The free energy (24) is the most general, universal and gauge i.e. frame rotation (11) invariant free energy that can be introduced to model the homopolymer in the limit where the characteristic length scales of spatial deformations around the minimum energy configuration become large in comparison to the distance (16) between consecutive monomers. It is essentially the naively discretized continuum Abelian Higgs Model with one complex scalar field, when expressed in the unitary gauge and the U(1) gauge transformation identified with the frame rotation (12); the term with parameter a is the Chern-Simons term which is commonly introduced in gauge theories to break parity. Thus, like the Abelian Higgs Model in general, it is unique in the sense of Kadanoff and Wilson. Implicit in (24) is the assumption that there are no abrupt wrenches and buckles along the polymer chain. Only small, gradual bends and twists are present in deviations around the energy minimum configuration, which is obtained by minimising the free energy (24). This is the limit of adiabatic deformations. The free energy (24) displays the discrete symmetry κi → −κi This symmetry may become spontaneously broken, by the ground state, in line with the conventional Abelian Higgs Model scenario; for the latter, see [30] in the present context. Moreover, the free energy (24) is not invariant under the local Z2 gauge symmetry (15): The free energy (24) is the leading non-trivial contribution to an expansion around a fixed background. To recuperate the Z2 symmetry one should interpret e.g. τ = 1 sin 2τ + O(τ 2 ) 2 etc. The free energy (24) can alternatively be interpreted as a deformation of the standard energy function of the discrete nonlinear Schr¨ odinger equation (DNLS) [31, 32]. The first two sums coincide the energy of the standard DNLS equation, in terms of the discretized Hasimoto variable [10]. The first (c) term in the third sum is the Proca mass that has a claim of gauge invariance; here it has been added as a ”regulator” for reasons explained in [13]. The second (a) term is the helicity, and the last (b) term is the conserved momentum. The last two terms break the Z2 parity symmetry, they are responsible for helicity of the homopolymer chain. In Table I the parameter values that are kept fixed during the simulations, have been listed. These parameter values have been chosen in conformity with those, that are commonly encountered in the case of proteins: For a protein, the torsion angles are much more flexible than TABLE I: The parameters in (24) that are kept fixed during our simulations. q 3.5 m 1.5 b 0 d 10−4 bond angles, and the value of m in Table I corresponds to an α-helical structure. The parameter a, which is not fixed, is of particular interest in the sequel. This is the parameter that breaks chirality. Note that the momentum of the DNLS hierarchy is not considered here i.e. b = 0, for simplicity. This term lacks a direct interpretation in the context of the Abelian Higgs Model. It turns out that the effects of this term are largely accounted for by the a dependent helicity, in any case. C. Pseudogap In [10] the following combination of the bond and torsion angles is considered ψi = σi eiϑi ≡ tan κi iϑi e 2 (25) where the phase is the anti-symmetrized combination ! i N X 1 X ϑi = τk − τk 2 k=1 k=i+1 In [10] it has been proposed to interpret (25) as the discretised version of the Hasimoto variable, in the Frenet frames. Accordingly ψi is the complex variable that converts (24) into a generalised version of the discrete nonlinear Schr¨odinger equation [31, 32]. When a symmetry becomes broken the effective potential (23) acquired a minimum where the modulus of the complex order parameter develops a non-vanishing expectation value, < |ψ| > 6= 0 (26) while in the unbroken, symmetric phase this quantity vanishes, at the level of the effective potential. The pseudogap state is a refinement of this phase structure [33– 37]. The system is in the pseudogap state when <ψ>= 0 (27) even though the modulus of ψ retains its non-vanishing symmetry breaking value (26), at the level of the effective potential. Such a situation occurs in the presence of a strong phase decoherence, < eiϑ > = 0 (28) 6 Consequently (27) can be observed even though the gap (26) is there. But the presence of a pseudogap can be difficult to detect, in terms of the effective potential alone, as the effective potential is commonly insensitive to the phase. A dynamical computation that engages fluctuations is needed. The pseudogap state is akin a symmetric phase precursor state in the broken symmetry phase. In particular the transition between the broken symmetry phase and the pseudogap state is commonly not a phase transition but a cross-over prelude to the fully symmetric state that enters when the lowest energy state of the effective potential is one where the modulus vanishes. It should be kept in mind, that (26)-(26) are commonly deduced by inspection of the effective potential. In the full theory, there are corrections due to fluctuations. In particular, in the full theory, at finite temperature, the modulus (26) never vanishes identically. The modulus is a positive definite quantity and thus, due to fluctuations, it always acquires a non-vanishing value in the full theory, as also shown in the simulations presented here. can annihilate each other, thus any pair of kinks can be removed by continuous deformations of the chain. A single kink can be translated, so that it becomes removed through the ends of the chain. However, it should be noted that on a discrete lattice the translation invariance is commonly broken, by the Peierls-Nabarro barrier [40–43]. Thus, in general it costs (thermal) energy to translate a kink along the chain. The present scenario is different from the one that appears in the case of kinks in folded proteins [38, 39]. There, the parameter values in (24) are different for different super-secondary structures i.e. helix-loop-helix motifs. A folded protein is described by a heteropolymer generalisation of (24), and the ground state is not a straight helix such as (31). A short analysis of a simple heteropolymer is presented in the sequel, in sub-Section IV E. E. D. Attractive interaction Zero temperature ground state The zero temperature ground state of the Hamiltonian (24) is a solution to the equations of motion, τi = a 1 2 2c dκi + 1 κi+1 = 2κi − κi−1 + 2q(κ2i − m2 )κi + cdτi2 κi (29) (30) Thus, with the present choice of parameters, the minimum energy ground state of (24) is a helix with q cd 2 τi ≈ ±m κi = ± m2 − 2q (31) a 1 a τi = 2c dm2 +1 ≈ 2c The parameter m has the fixed value given in Table I throughout. As a consequence the pitch of the ground state helix is controlled by the ratio a/c. In the sequel the effect of helicity on the phase structure is investigated by varying the ratio a/c within the range a/c ∈ [0, 4π]. For some parameter values, there can be a conflict between the values of (κi , τi ) given by (31) and the steric constraint (17). In that case the lowest energy ground state configuration is obtained by minimising (24), subject to the constraint (17) The range of the bond angle has been extended to negative values, by the Z2 symmetry (15). The two ground states κi = ±m have the same energy. In addition of these two ground states, there can also be local minima of (24) that have the profile of a kink [38, 39] i.e. a domain wall that interpolates between the two ground states κi = ±m. The energy of a kink is higher than the energy of the ground state helix (31). Two kinks In the presence of an ambient temperature the bond and torsion angles fluctuate. This gives rise to a nonvanishing amplitude for finite energy barrier penetration. Thus, thermal fluctuations eventually remove any kink there is, and the system reaches a thermal equilibrium which minimises the Helmhotz free energy. At very low temperature values, this is a configuration which is very close to the ground state of (24), (17). At very low temperatures, the chain is then essentially a straight, linear helical rod. On the other hand, in the limit of very high temperatures the constraint (17) dominates. The chain enters the self-avoiding random walk phase, characterised by the mean field compactness index value ν = 3/5. Between the low temperature rod and high temperature SARW limits there is a priori a transition region. Recall that the conventional θ-regime is located between the collapsed phase and the self-avoiding random walk phase [5]. Thus this transition regime between a rod and SARW, if it indeed exists, may be called the η-regime. Among the goals of the present investigation is to search for such an η-regime, and the relation between the θ- and η-regimes and a pseudogap state. It should be noted that the constraint (17) is purely repulsive. There are no attractive interactions present in the model. For a collapsed phase, an attractive interaction between the monomers must be introduced [5]. In particular, the θ-regime is a consequence of a long distance balance between attractive interactions and repulsive interactions, to the effect that no long-range interactions prevail [5]. Thus, for a collapsed phase, the model (17), (24) is amended with an attractive interaction. Accordingly, (17) is generalised into an interaction potential with the 7 following profile, U (r) = +∞ 0<r<δ U (r) = U0 {tanh(r − R0 ) − 1} δ < r < +∞ (32) Here δ is the radius of the self-avoiding condition (17). Hence for r < 0 the steric constraint (17) persists. But for r > δ there is a short range attractive interaction with strength determined by the parameter U0 . In the sequel this parameter will be varied, jointly with the ratio a/c that characterises helicity. The parameter R0 which determines the range of the attractive interaction, has the following value R0 = 5.0 (˚ A) throughout. It is reminded, that in all-atom molecular dynamics simulations any long range interaction between atoms is commonly cut off, sharply, beyond distances around 10 ˚ Angstr¨ om or so. The attractive interaction has a physical interpretation, in terms of ”hydrophobic” forces: In the case of e.g. a protein chain under physiological conditions, there is an effective attractive interaction between those amino acids which are considered ”hydrophobic” [5–7]. Thus in the presence of the attractive interaction (32) the energy function (24) models a chain made of ”hydrophobic” residues, with ”hydrophobicity” that depends on the value of U0 . It is noted that qualitatively, the present results have been found to be quite insensitive to the details of the profile of the potential U (r). III. SIMULATION DETAILS Extensive simulations have been performed, to test three different algorithms [44]. The steric constraint (17) has been used in lieu of the more general (32). The aim is to identify a Monte Carlo algorithm that has the fastest rate of convergence towards a thermal equilibrium state. The three Markovian algorithms that have been tested are 1) Heat Bath algorithm 2) Metropolis algorithm 3) Mixed algorithm The curvature and torsion angles are updated according to a probability distribution, that satisfies the detailed balance condition of a canonical ensemble is obtained in the limit of an infinite number of updates. Each update consists of a ”walk” through the entire chain with a provisional revision of each value (κi , τi ). The three algorithms differ from each other only in the manner how the new values , τinew ) are generated. (κnew i 1. Heat Bath algorithm In the Heat Bath algorithm, new values (κnew , τinew ) i are generated randomly, according to probability distributions 1 P (κnew ) = exp {−βHi,κ (κnew )} (33) i i Zi,κ and P (τinew ) = 1 exp {−βHi,τ (τinew )} Zi,τ (34) Here Hi,κ and Hi,τ are the sum of all those terms in the Hamiltonian (24) that contain κi resp. τi with the given index i and the Zi,κ and Zi,τ are normalisation factors. The updated values of κnew and τinew do not depend on i the previous values of κi and τi . The probability density for κnew has the form i P (κi ) ∼ exp{−c1 κ4i − c2 κ2i − c3 κi } (35) where c1 = βq c c2 = β(2 − 2qm2 + dτi2 − abτi ) 2 c3 = β(−2(κi+1 + κi−1 )) (36) Thus (35) is non-Gaußian. On the other hand, the probability density P (τinew ) has the Gaußian profile c P (τi ) ∼ exp{−β( [dκ2i + 1]τi2 − aτi )} 2 (37) Rejection sampling has been used to generate random numbers according to these probability distribution: After generating κnew and τinew the self-avoidance condii tion (17) is checked, and the update is rejected when the condition is violated. 2. Metropolis algorithm Here β is the inverse Monte Carlo temperature. The equilibrium distribution New values of κ and τ are generated according to Gaußian probability distributions, which is centered at the old values. The dispersion of each Gaußian can be adjusted, to enhance the convergence of the algorithm. The new values of κ and τ are accepted or rejected, in the same manner as in the conventional Metropolis algorithm. For example, in the case of τi the probability of acceptance of a new value is exp{−βH(κ, τ )} P (τi ) = min{1, exp(−β∆H)} P ({κnew , τnew }, {κold , τold }) exp(−βH({κold , τold })) = P ({κold , τold }, {κnew , τnew }) exp(−βH({κnew , τnew })) (38) 8 where ∆H is the difference of the energy between the new and the old configurations. In addition, the selfavoidance condition (17) is also verified at each step. Mixed algorithm The values of κnew are generated in the same mani ner as in the Heat Bath algorithm, while for τinew the Metropolis algorithm is used. The convergence of the algorithm can be adjusted, by changing the dispersion of the Gaußian distribution in the τ update. 80 Metropolis Mixed, ∆τ = 0.02 Mixed, ∆τ = 0.20 Mixed, ∆τ = 0.50 Mixed, ∆τ = 1.00 Mixed, ∆τ = 2.00 Heat Bath 70 60 Rgyr 3. 50 40 30 A. Comparison of algorithms 20 10 A series of comparative simulations have been performed, using the three algorithms. The following parameter values have been used in the hamiltonian (24), (39) These values have been chosen to reproduce a monotonous α-helical structure of a protein, as the lowest energy conformation; the choice is not unique. The simulations have been performed with varying chain lengths, from N=100 to N=900. In each case, simulated annealing has been used and the initial configuration is always a linear straight rod with κi = τi = 0 -10 -8 -6 -4 log10T -2 0 2 4 FIG. 1: Dependence of Rgyr on temperature for identical homopolymer chains simulated with different algorithms. Here ∆τ is dispersion of Gaußian distribution for generation of τnew in mixed algorithm. In conventional Metropolis algorithm the dispersions are ∆κ = ∆τ = 0.01 45 n = 100 n = 200 n = 300 n = 400 40 35 Rgyr a = −1.0 × 10−4 c = 1.0 × 10−4 d = −5.8 × 10−4 -12 30 25 20 The initial configuration is first heated to very high temperature values (up to Tmax = 100), where the structure is fully thermally randomised. This is followed by a slow cooling period, to the target temperature. The cooling takes place with small temperature steps, with each step equal to ∆T = 0.05 in logarithmic scale. After each step, 103 − 5 × 104 Monte Carlo updates are performed along the whole chain, to ensure that it becomes thermalised to the ambient temperature; the longer the chain, the longer the thermalisation. The final phase diagrams have been calculated using 5000 Monte Carlo updates per one step of cooling. In Figure 1 the value of the radius of gyration (1) is compared, as a function of Monte Carlo temperature for a chain with N=100 monomers. The results obtained using the Metropolis algorithm are found to be very different from those obtained using the Heat Bath and Mixed algorithms; for small dispersion ∆τ the Mixed algorithm coincides with the Metropolis algorithm. But when ∆τ increases, the Mixed algorithm approaches the Heat Bath algorithm as shown in the Figure. There is no a priori reason why the results for the three algorithms should be different: The stationary distribution is the same, in each of the three algorithms. But apparently one of the three converges very slowly 15 10 10 100 1000 10000 Termalization length 100000 FIG. 2: Dependence of Rgyr on the thermalisation length in Metropolis algorithm. The number of updates, per one step in simulated annealing process, is shown along the horizontal axis. towards the equilibrium distribution. The Markovian processes have been modelled by simulations, which are finite in time. Thus, in order to conclude which of the algorithms is slow to converge, the dependence of the result on the length of simulation needs to be analysed. In the case of the Metropolis algorithm, the results are shown in the Figure 2, for T = 10−13 which is the lowest Monte Carlo temperature value that has been used in the present simulations. As shown in the Figure, Rgyr continues to increase with increasing length of thermalisation. Thus, it is concluded that the results obtained with the Metropolis algorithm apparently approach those of the Heat Bath algorithm very slowly as a function of the simulation time. Accordingly, it is concluded that the 9 200 35 log T = 2 log T = -2 log T = -10 160 n = 1000 n = 2000 n = 5000 n = 10000 30 25 Rgyr Rgyr 120 80 20 15 ν2 = 0.593 ± 0.022 ν-2 = 0.622 ± 0.014 ν-10 = 0.917 ± 0.073 40 0 200 400 600 800 10 5 1000 0 0.5 Metropolis algorithm is the least optimal for the present purposes. Either the Heat Bath algorithm, or alternatively the Mixed algorithm with sufficiently large ∆τ , should be used, to try and describe the thermal equilibrium configurations. In addition, the compactness index ν i.e. the inverse of the Hausdorff dimension has been inspected using the three different algorithms. The results for the Heat Bath algorithm are shown in Figure 3. Between T = 100 and T = 1000 only very small variations are detected in the value of the radius of gyration Rgyr ; the value of ν is essentially temperature independent, and apparently corresponds to the SARW phase. At very low temperatures a transition to the rigid rod state, with ν ≈ 1, is observed. This result is in line with the general arguments that have been presented in sub-section II E, on the expected phase structure of the homopolymer model (24), (17). IV. COLLAPSED PHASE A. The algorithm Extensive Monte Carlo simulations have been performed with the model (24) when subject to the attractive interaction (32). On the basis of the test simulations, described in sub-section III A, the following algorithm is proposed: The probability distribution is P = X 1 exp{−β[H + U (~ri − ~rj )]} Z i<j Here H is the Hamiltonian (24) and U is the potential (32). The Metropolis algorithm is used for acceptance, but with a proposal distribution that coincides with the Heat Bath algorithm. That is, new values of κi and τi are generated using the distributions (33) and (34). The ensuing homopolymer configuration is then accepted, provided it satisfies both the self-avoidance condition (17) 1.5 2 2.5 log10T N FIG. 3: Compactness index ν at different temperatures for a homopolymer chain. Simulation performed with Heat Bath algorithm 1 FIG. 4: Dependence of Rgyr on temperature T for a chain with 100 monomers, in the truncated model. Here n is the number of Monte Carlo updates per one simulated annealing step in temperature. and the Metropolis accept-reject condition that utilises the residual energy X EU = U (~ri − ~rj ) i<j The acceptance criterion is exp(−β∆EU ) > λ where λ is a random number which is uniformly distributed between 0 and 1, and ∆EU is the change in EU under the update of κi and τi . B. Truncated model The algorithm has been tested in a truncated model, where the Hamiltonian H is removed and only the contribution from the potential U is retained. The radius of gyration Rgyr has been evaluated, as a function of the Monte Carlo temperature. The results are shown in Figure 4 for several thermalization lengths. The value U0 = 15 is used in (32) for all the runs; the exact numerical value is not relevant in these simulations, it can be adjusted by adjusting the scale of MC temperature T A smooth transition is observed in the radius of gyration, from larger values at high temperatures to smaller values at low temperatures. The Figure 5 shows how Rgyr depends on the length of polymer chain N , and fitted to the leading order contribution in (2). It is found that Rgyr decreases as the thermalisation length increases. The value of the compactness index ν also decreases, and converges towards the mean field value 1/3 of the collapsed phase, as shown in Figure 6. For high temperatures a ν ≈ 3/5 is found, corresponding to the SARW phase. It is concluded that the truncated model displays two phases: At high temperatures there is the self-avoiding random walk phase (SARW) and at low temperatures there is the collapsed phase. 10 20 18 16 14 Rgyr mer enters a regime of decreasing Rgyr . Finally, there is the low temperature regime where the radius of gyration Rgyr has a small value; see Figure 9. The compactness index ν shows that when the thermalisation increases the value of ν converges towards ν ≈ 1/3. This is shown in Figure 10. n = 1000 n = 2000 n = 5000 n = 10000 12 10 ν1000 = 0.431 ± 0.012 ν2000 = 0.3997 ± 0.0079 ν5000 = 0.3668 ± 0.0011 ν10000 = 0.3679 ± 0.0036 8 6 4 100 200 300 400 500 2. 600 N FIG. 5: Dependence of Rgyr on the length of polymer chain N for different thermalisation length with T = 1, in the truncated model. According to (31) the classical ground state profile of τi remains intact when the parameters a and c are changed in such a manner that the ratio a/c is constant. To study the effect of such a change in a and c at finite temperature, in particular how it conspires with the parameter U0 , simulations have first been performed with a = c = 10−2 > U0 = 10−4 0.46 0.44 ν 0.42 0.4 0.38 0.36 2000 4000 6000 8000 Length of thermalization 10000 FIG. 6: Dependence of compactness index ν on thermalisation length when T = 1, in the truncated model. C. Parameters a and c Full model and parameter effects The variant of the Heat Bath algorithm described in sub-section IV A has been used to study the effect of various parameters to the phase structure of the full homopolymer model with attractive interaction (24), (32). These results are summarised in Figures 7-11. 1. Parameter U0 Figures 7-9 describe the properties of the radius of gyration for three representative values of the strength parameter U0 in (32), 10−4 10−2 U0 = 10−1 For the other parameters the values (39) are used. For each value of U0 , three different characteristic regimes are observed. In the high temperature limit the homopolymer is found in the SARW phase; this is confirmed in Figure 8. When the temperature decreases, the homopoly- (40) with the values (39) for the remaining parameters. Thus, unlike in the previous simulations now the characteristic scale of the attractive interaction is smaller than that of the torsion angle dependent terms in the Hamiltonian. The results for the radius of gyration are presented in the Figure 11. At high temperatures the chain is again in the SARW phase. Then, as temperature decreases, there is a transition to a regime akin the intermediate regime shown in Figure 7. Finally, there is a low temperature regime where the chain fluctuates around the classical solution (31). The scale of transition to the low temperature regime is controlled by parameters a and c in Hamiltonian. It is concluded that when the scale U0 of the shortrange attractive interaction is smaller than the scale of the parameters a and c, the low temperature limit is described by helical structures. D. Analysis of different phase regimes The bond and torsion angles form the complete set of local order parameters to probe the phase structure (3), in the case of the present homopolymer model. These order parameters have the following characteristics, in the different regimes that have been analysed in Figures 7-11; the results are summarised in Figures 12-15. In the very high temperature SARW phase, both the bond angle and the torsion angle are subject to large fluctuations; the simulation results are shown in Figure 12. For the bond angles, the values are distributed in the range 0 ≤ <|κ|> ≤ κmax ∼ 2.2 where the upper limit is due to the steric constraint (17). As temperature increases, the values of κ become increasingly evenly distributed over this range so that in the T → ∞ limit the distribution is fully uniform; the Figure 12 shows the bond and torsion angle distribution at a 11 24 U = 0.0001 U = 0.01 U = 0.1 22 20 50 45 40 16 Rgyr Rgyr 18 14 35 12 30 10 25 8 6 -8 -6 -4 -2 0 20 2 -8 -6 -4 -2 0 2 log10T log10T FIG. 7: Dependence of Rgyr on temperature for a chain with 50 monomers. Thermalisation length is n = 104 updates per one step of simulated annealing FIG. 11: Dependence of Rgyr on temperature. The chain has 100 monomers, and parameters (a, c, U0 ) are given in (40). 90 70 Rgyr 60 50 40 30 20 ν = 0.597 ± 0.019 10 0 100 200 300 N 400 500 600 n = 1000 n = 2000 n = 5000 n = 10000 n = 50000 24 Rgyr 20 16 8 4 100 200 300 N 400 500 8 4 0 -2.5 -1.5 -0.5 0.5 1.5 12 8 4 0 2.5 -7 -4.2 -1.4 1.4 4.2 7 τ FIG. 12: Distribution of (κi , τi ) in the high temperature SARW phase. generic but high temperature value. It is apparent from this Figure that both angles are disordered. It is concluded that the SARW phase is a disordered phase. A priori, the intermediate transition regime should be the canonical θ-regime; the transition is found to take place between temperature values within the range 10−3 < T < 101 ν1000 = 0.4719 ± 0.0045 ν2000 = 0.442 ± 0.012 ν5000 = 0.4327 ± 0.0058 ν10000 = 0.434 ± 0.015 ν20000 = 0.3932 ± 0.0013 12 12 16 κ FIG. 8: Dependence of Rgyr on the length of chain for high temperature T = 101.7 and with U0 = 0.01. 28 16 Probability distribution Probability distribution 80 600 as can be seen in Figures 7 and 11. In this putative θregime the values of the bond angle are found to become ordered. This is shown in Figure (13): The values of κi are thermally fluctuating at around κ ≈ 1.5. The values of the torsion angle remain largely disordered. But the 0.48 0.46 ν 0.44 0.42 0.5 κ 1.5 2.5 16 12 8 4 0 -7 -4.2 -1.4 1.4 4.2 τ FIG. 13: Distribution of (κi , τi ) in the intermediate transition regime. 0.4 0.38 48 40 32 24 16 8 0 -2.5 -1.5 -0.5 Probability distribution Probability distribution FIG. 9: Dependence of Rgyr on the length of chain for low temperature T = 10−6 and various thermalisation lengths, and with U0 = 0.01. 1000 10000 Length of thermalization 100000 FIG. 10: Dependence of compactness index ν on thermalisation length. Temperature is equal to T = 10−6 and U0 = 0.01. 7 12 60 40 20 0 -2.5 -1.5 -0.5 0.5 1.5 25 40 15 35 10 5 30 25 0 2.5 n = 1000 n = 2000 n = 5000 n = 10000 45 20 Rgyr Probability distribution Probability distribution 50 80 -7 -4.2 -1.4 κ 1.4 4.2 7 20 τ 15 FIG. 14: Distribution of (κi , τi ) in the collapsed phase. 10 -8 -6 -4 -2 0 2 80 Probability distribution Probability distribution log10T 60 40 20 0 -2.5 -1.5 -0.5 0.5 κ 1.5 2.5 80 FIG. 16: Dependence of Rgyr on temperature for a heteropolymer chain with 150 monomers Thermalization length n is number of updates per one step of simulated annealing 60 40 20 16 0 -2 -1 0 1 2 15 τ 14 FIG. 15: Distribution of (κi , τi ) in the helical rod-like phase. Rgyr 13 emergence of a structure which is indicative of a double well effective potential is observed to be present in the τ -distribution; see Figure 13. The intermediate region is identified as a pseudogap state, as described in sub-Section II C. Finally, there are two different low temperature phases: The collapsed phase shown in Figure 7 where the attractive short-distance interaction dominates and the helical rod-like phase shown in Figure 11 where the attractive short-distance interaction becomes weak. The distribution of bond and torsion angles in the collapsed phase are displayed in Figure 14. The bond angle is highly ordered around the classical value (31) but the torsion angle remains disordered. However, there is an apparent spontaneous symmetry breaking that has taken place; the double well structure seen in the τ distribution of Figure 13 has been removed, in a way that resembles the familiar spontaneous symmetry breaking in a Z2 symmetric potential well. In the helical rod-like phase, both the bond and torsion angles become peaked around the classical values (31). The configurations are akin straight helical rods; a little like e.g. collagen when biologically active. The κ distribution reflects the discrete Z2 gauge symmetry. But the Z2 symmetry in the τ distribution observed in Figure 13 is fully broken. The transition between the collapsed phase and the helical rod-like phase entails a transitition where the torsion angles become ordered. Due to the very low temperature values involved, the fluctuations are strongly suppressed and a simulation becomes tedious. It is conjectured that, when the temperature is kept in the low temperature regime, an initial helical rod-like structure but with parameter values corresponding to the collapsed phase, is in 12 11 10 ν = 0.269 ± 0.039 9 8 50 100 150 200 250 N 300 350 400 FIG. 17: Dependence of Rgyr on the length of polymer chain for low temperature T = 10−7.5 and various thermalisation lengths. a glassy phase. Vice versa, an initial collapsed configuration with parameter values in the helical rod-like phase, will eventually become subject to cold denaturation. E. heteropolymer The effect of perturbations that break the homogeneity of homopolymer have been investigated as follows: A collapse has been found to take place when the parameter U0 that characterises the strength of self-interaction is larger than the parameters a and c that characterise the torsion angle dependent terms. Thus, a short segment is introduced along the chain, where U0 is less than a and c. Accordingly, a simulation is performed where a heteropolymer is constructed so that for a short sub-chain of 12 monomers, the values of the parameters a and c is increased from a = c = 10−6 to a = c = 10−2 , while U0 = 10−4 along the entire chain. Results of simulations are presented in Figures 16 and 17, for a chain with 150 monomers. The dependence of Rgyr on temperature is shown in figure 16; the phase diagram is very similar to 7. For example, at low temperatures it is found that 13 -5.4 -5.7 -6 -2 log10U0 -3 -4 -5 -6 FIG. 18: (left) Homopolymer does not display helices in the ground state. (right) Addition of a subchain that breaks homogeneity, gives rise to a helix in the subchain region. 4 3 a/(πc) 2 1 0 -7 -5 -3 -1 1 log10T FIG. 19: The phase diagram on the (U0 , T, a) space. the compactness index is close to the mean field value ν = 1/3 of the collapsed phase. However, the geometry of a configuration in the collapsed phase is different: As shown in Figure 18 a helical structure appears only in that sub-chain where the parameter values a and c have been increased. V. PHASE DIAGRAM The homopolymer phase is found to depend on three relevant scales. – There is the extrinsic temperature scale where the values of κi become ordered. This scale can also be controlled intrinsically, by the parameter q in (24), but the details have not been addressed here. – There is temperature scale where the values of τi become ordered. This scale can be controlled intrinsically, by the parameter ratio a/c in (24). – Finally, the effects of the scale U0 for the short-range attractive interactions have been investigated. This parameter determines an intrinsic scale that controls the transition temperature alternatively to the collapsed phase, or to the helical rod-like phase. Thus, by changing the relations between the three scales the phase diagram of the homopolymer can be identified; the phase diagram is constructed here in terms of (T, a, U0 ). All the remaining parameters are fixed, and given by the values in Table I. The Figure 19 shows the three-dimensional phase diagram in the (T, a, U0 ) space. The Figures 20-23 show various cross-sections, taken at selected values of U0 i.e. these Figures show the phase diagram in Figure 19 on the (a, T ) plane, with different values of U0 . Figure 20 shows the phase diagram, for a quite large value of U0 (strong coupling) At high temperatures there is the SARW phase. When temperature decreases there is the pseudogap state, that becomes the collapsed phase at low temperatures. Between the pseudogap state and the low temperature collapsed phase there is a θ-regime, or rather a θ-point as it is observed only over a very narrow temperature range. • Figure 21 displays the phase diagram, when the value of U0 is lowered but still relatively large (intermediate but not weak coupling). At high temperatures there is again the SARW phase, followed by the pseudoogap state as the temperature decreases. At low temperatures the pseudogap state becomes converted either to the collapsed phase or to the straight rod phase, depending on the value of the helicity parameter a. In addition, there is a range of values of a, when a intermediate similar to the θ-regime is observed between the pseudogap state and the straight rod phase. This is the η-regime. It is notable, that there is a possibility of a 4-critical point involving the pseudogap state, the η-regime, and the collapsed and straight rod phases. • Figure 22 shows the phase diagram, as the value of U0 becomes further decreased (intermediate but not strong coupling). The collapsed state has entire disappeared, and replaced by the straight rod phase at very low temperatures. The η-regime displays a periodic structure in the parameter a. There appears to be a tri-critical point involving the pseudogap state, the η-regime and the straight rod phase. • Finally, in Figure 23 the weak coupling U0 phase diagram is displayed. The overall topology of the phase diagram is similar to the one in Figure 20, but with the straight rod phase as the low temperature phase instead of the collapsed phase. The Figure 24 shows the phase diagram on the (U0 , T ) plane, with helicity fixed. It is notable that there might be a 5-criticality involving the η- and θ-regimes, the pseudogap state and the collapsed and straight rod phases. However, the detailed investigation of this region of the phase diagram is beyond the capacity of the computer power which is presently available to us. In summary, it is remarkable how the phase diagram 14 2 SARW -3 pseudogap -2 2.5 -2.5 log10U = -3 0 2 θ-regime -4 collapsed -6 log10U0 log10T -3.5 e re θ- -4 -4.5 1.5 gim collapsed 1 pseudogap 0.5 SARW -5 -8 0 0.5 1 1.5 2 a/(πc) 2.5 3 3.5 0 straight -5.5 rod -6 -6 η-regime -5 -4 FIG. 20: A cross-section of the phase diagram in 14 at U0 = 10−3 . -0.5 -3 1 2 -1 log10T FIG. 24: The phase diagram on the (U0 , T ) plane, with a = 10−4 . 2 SARW log10U = -4.8 0 pseudogap log10T -2 collapsed -4 η-regime straight rod -6 -8 0 0.5 1 1.5 2 a/(πc) 2.5 3 3.5 FIG. 21: A cross-section of the phase diagram in 14 at U0 = 10−4.8 .. 2 SARW log10U = -5.1 0 pseudogap log10T -2 η-regime -4 -6 straight rod -8 0 0.5 1 1.5 2 a/(πc) 2.5 3 3.5 (19)-(23) is periodic in the helicity parameter a (when c is fixed). Moreover, there is a rapid transition between the helical phase where the compactness index ν = 1 and the collapsed phase where ν ≈ 1/3 at low temperatures, as shown in Figure (21). There is a pseudogap state that appears as a transition regime, akin the conventional θ-regime, between the collapsed and the SARW phases. It is conjectured that in the thermodynamical limit where the number of monomers N → ∞ the value of the compactness index in this pseudogap state becomes ν = 1/2. For the transition regime between SARW and helical phases shown in Figure (22) there is a pseudogap state that should essentially coincide in its properties with the θ-regime pseudogap state. But due to computational limitations the present analysis is not sufficient to confirm this. There is an apparent 4-critical, even 5-critical point as shown in Figures 21 and 24, but the detailed analysis of this region in the phase diagram needs to be performed using more extensive simulations, which is postponed to a future project. However, it is observed that the potential presence of a 4-critical point in a theoretical context very similar to the present one (i.e. Abelian Higgs model) has been reported in [26] FIG. 22: A cross-section of the phase diagram in 14 at U0 = 10−5.1 . VI. 2 SARW log10U = -6 0 pseudogap log10T -2 -4 η-regime -6 straight rod -8 0 0.5 1 1.5 2 a/(πc) 2.5 3 3.5 FIG. 23: A cross-section of the phase diagram in 14 at U0 = 10−6 SUMMARY The phase structure of chiral homopolymers under varying ambient temperature values have been investigated theoretically, in terms of a universal infrared limit energy function in combination with steric constraints (self-avoidance) and a short-range attractive interaction between residues. As such, the model should provide a realistic description even in the case of heteropolymers that display an approximatively repeating monomer pattern, provided the scale of the repeat can be considered small in comparison to the chain length. A biologically important example is given by collagen, the most prevalent protein in our body, in which case the short range attractive interaction models weak hydrophobicity of the amino acids. 15 It is found that the phase diagram displays a high level of complexity, in terms of the parameters that control the helicity, and the strength of the attractive interaction. In particular, the low energy phase is either like a linear one dimensional straight rod, or a space filling collapsed configuration. At intermediate temperatures, there is a state which can be identified as an example of the pseudogap state, and there are also intermediates that are more like the conventional θ-regime. It is possible that these regimes merge, in the thermodynamical limit, at least for some range of parameter values. However, the possibility of the existence of 4-critical, even 5-critical points in the phase diagram is also proposed, but can not be confirmed with presently available computer power. The extension of the present approach to investigate properties of proteins and other heteropolymers remains a challenge to future research. [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] A.J. Heeger, Rev. Mod. Phys. 73 681 (2001) A.J. Niemi,G.W. Semenoff, Phys. Repts. 135 99 (1986) J. Hiraki, Journ. Antib. Antif. Ag. 23 349 (1995) C. Li, D.F. Yu, A. Newman, F. Cabral, C. Stephens, N.R. Hunter, L. Milas, S. Wallace Cancer Res. 58 2404 (1998) P.G. DeGennes, Scaling Concepts in Polymer Physics (Cornell University Press, Ithaca, 1979) A.Yu. Grosberg, A.R. Khokhlov Statistical Physics of Macromolecules (AIP Series in Polymers and Complex Materials, Woodbury, 1994) L. Sch¨ afer, Excluded Volume Effects in Polymer Solutions, as Explained by the Renormalization Group (Springer Verlag, Berlin, 1999) A.J. Niemi, Phys. Rev. D67 106004 (2003) U.H. Danielsson M. Lundgren, A.J. Niemi, Phys. Rev. E82 021910 (2010) S. Hu, Y. Jiang, A.J. Niemi, Phys. Rev. D87 105011 (2013) T. Ioannidou, Y. Jiang, A.J. Niemi, Phys. Rev. D90 025012 (2014) A.J. Niemi, Theor. Math. Phys. 181 1235 (2014) A.J. Niemi, arXiv:1412.8321 [cond-mat.soft] L.P. Kadanoff, Physics 2 263 (1966) K.G. Wilson, Phys. Rev. B4 3174 (1971) K.G. Wilson, J. Kogut, Phys. Repts. 12 75 (1974) B. Widom, J. Chem. Phys. 43 3892 (1965) M.E. Fisher, Rev. Mod. Phys. 46 597 (1974) P.G. De Gennes, Phys. Lett. 38A 339 (1972) J.C. LeGuillou, J. Zinn-Justin, Phys. Rev. B21 3976 (1980) B. Li, N. Madras, A. Sokal, Journ. Stat. Phys. 80 661 (1995) A. Krokhotin, A. Liwo, A.J. Niemi, H.A. Scheraga Journ. Chem. Phys. 137 035101 (2012) M.L. Huggins, Journ. Chem. Phys. 9 440 (1941) P.J. Flory, Journ. Chem. Phys. 9 660 (1941) VII. ACKNOWLEDGEMENTS AJN acknowledges support from Region Centre Recherche d0 Initiative Academique grant, Vetenskapsr˚ adet, Carl Trygger’s Stiftelse f¨or vetenskaplig forskning, and Qian Ren Grant at BIT. The work of MU was supported by the DFG grant SFB/TR-55 and by Grant RFBR-1402-01261-a. Numerical calculations were performed at the ITEP computer systems “Graphyn” and “Stakan” and at the supercomputer center of Moscow State University. AJN thanks A. Sieradzan and Nevena Ilieva for discussions. We dedicate our research to the memory of M. Polikarpov. A close colleague, mentor and teacher without whom our collaboration would not be. Thank you, Misha. [25] P.J. Flory, Principles of Polymer Chemistry (Cornell University Press, Ithaca, 1953) [26] M. Bock, M.N. Chernodub, E.-M. Ilgenfritz, A. Schiller Phys.Rev. B76 (2007) 184502 [27] S. Hu, M. Lundgren, A.J. Niemi, Phys. Rev. E83 061908 (2011) [28] K. Hinsen, S. Hu, G.R. Kneller, A.J. Niemi Journ. Chem. Phys. 139 124115 (2013) [29] S. Coleman, E. Weinberg Phys. Rev. D7 1888 (1973) [30] M.N. Chernodub, L.D. Faddeev, A.J. Niemi, JHEP 12 014 (2008) [31] L.D. Faddeev, L.A. Takhtajan, Hamiltonian methods in the theory of solitons (Springer Verlag, Berlin, 1987) [32] M.J. Ablowitz, B. Prinari, A.D. Trubatch, Discrete and Continuous Nonlinear Schr¨ odinger Systems (London Mat. Soc. Lect. Note Series 302, London, 2003) [33] V.J. Emery, S.A. Kivelson, Nature 374 434 (1995) [34] J. Corson, R. Mallozzi, J. Orenstein, J.N. Eckstein, I. Bozovic, Nature 398 221 (1999) [35] N. Mannella, et.al. Nature 438 474 (2005) [36] H. Kleinert, E. Babaev, Phys. Lett. B438 311 (1998) [37] K. Zarembo, JETP Letters 75 59 (2002) [38] M. Chernodub, S. Hu, A.J. Niemi, Phys. Rev. E82 011916 (2010) [39] N. Molkenthin, S. Hu, A.J. Niemi, Phys. Rev. Lett. 106 078102 (2011) [40] R. Peierls, Proc. Phys. Soc. 52 34 (1940) [41] F. Nabarro, Proc. Phys. Soc. 59 256 (1947) [42] F. Nabarro, Mat. Sci. Eng. A234 67 (1997) [43] A.K. Sieradzan, A.J. Niemi, X.Peng, Phys. Rev. E90 062717 (2014) [44] B.A. Berg, Markov Chain Monte Carlo Simulations and Their Statistical Analysis (World Scientific, Singapore, 2004)
© Copyright 2024