Physical mechanism for biopolymers to aggregate and maintain in non-equilibrium states

Many human or animal diseases are related to aggregation of proteins. A viable biological organism should maintain in non-equilibrium states. How protein aggregate and why biological organisms can maintain in non-equilibrium states are not well understood. As a first step to understand such complex systems problems, we consider simple model systems containing polymer chains and solvent particles. The strength of the spring to connect two neighboring monomers in a polymer chain is controlled by a parameter s with s → ∞ for rigid-bond. The strengths of bending and torsion angle dependent interactions are controlled by a parameter s A with s A → −∞ corresponding to no bending and torsion angle dependent interactions. We find that for very small s A, polymer chains tend to aggregate spontaneously and the trend is independent of the strength of spring. For strong springs, the speed distribution of monomers in the parallel (along the direction of the spring to connect two neighboring monomers) and perpendicular directions have different effective temperatures and such systems are in non-equilibrium states.

parts of equal probable microscopic states. The non equilibrium processes of biological relevance often, however, probes only very limited parts of those states, which involve collective properties of limited number of molecules under spatial and temporal constrictions. The formation of amyloid-like fibrils from peptides or protein molecules is one of such processes, which is currently in focusing for biomedical researches 11,22,[28][29][30][31][32] and for potential material applications 33 . Over the aggregation process, the accumulated domains are composed of lined-up segments of β-strands 34 where the stacking can be enhanced by external confining 35 . The geometric packing, with the assistance of hydrogen-bonding, therefore, plays an important role in the process. Since the efficient stacking of individual protein molecule is essentially guided by the backbone-connectivity 36,37 , the focus on the latter property and on the combined effects with external confinement may provide the most relevant factors for the processes shared by a wide range of different peptides. Such issues can always be addressed efficiently by studying simplified models [38][39][40] .
One important issue is on the nature of the process as a nucleation mechanism with its origin from thermodynamic instability. The mechanism is known to be responsible for the ordering of structured chains of homopolymer 41 and the gelation for stacking of spherical molecules 42,43 . Some recent experimental results for protein aggregation 44,45 are compatible with such an aspect. To focus on those features 35,44 shared by the solidification process of simple molecules, in either bulk of small volume 46 or in confined space 47 , and meanwhile to underscore the effects of backbone connectivity in a systematic manner, we have initiated studies by molecular dynamics simulations 48,49 , on the aggregation problem of dilute polymer chains with effective soft confinement imposed by the finite space and the boundary conditions in simulation. While the main bodies of the clusters formed in these simulations are in bundled structures, which are different from the ordered layers obtained in processes of crystallization in polymer melt 50 , the dynamics of ordering of the dilute chains can still be quantified effectively 48 . It is also shown 48, 49 that we can manipulate the ruggedness along the backbones to prepare systems of polymer chains which intrinsically do not undergo processes of aggregation. For these systems, enhanced stiffness may render the chains stay at non-equilibrium 'quasi-steady' situations, which show unusual thermal equilibration properties in presence of contact with background fluid and have non Maxwell-Boltzmann velocity distributions for the monomers 51,52 . One question of particular interest is then whether this kind of stability implicitly maintained by the backbone connectivity is also present in the middle of aggregation processes occurring for those systems which are subject to thermodynamic instability.
The model systems under our consideration are in fixed volume, with their steady non equilibrium phases within the thermodynamically unstable regime (Fig. 1) determined by the strength parameter K A of angle potential along the chains in powers of ten, = K 10 s A A . The collection of chains will have a bundled structure after going through the process from stage I to stage IV ( Fig. 1(d)), reaching steady phases in the lower temperature regime of Fig. 1(e), if K A of the system is sufficiently small (in yellow colored region). We have chosen a specific minimum for the bending angle potential (see Methods) to introduce local hindrance into individual chains, which prevents the backbones from forming lined-up segments and the ordering is prohibited as soon as the strength K A is increased to push the system entering the grey colored region in Fig. 1(e). Under the circumstance the stage IV can never be reached and the system is left with configurations of partially ordered local conformations. The occurrence of ordering is insensitive to the bonding strength k n.n. and is enhanced but not determined by the spatial confinement 48,49 .
In previous studies 51, 52 , we have revealed in simulations, the primitive form of backbone induced stability for systems having s A = −1, which undergo no clustering at all temperatures ( Fig. 1(e)). The simulations have been carried out without imposing external temperature controls and the numerical noises in the simulations are treated as unavoidable sources of heat, causing the systems out of equilibrium, with temperatures varied in time (along the red dashed line on the right of Fig. 1(e)). The backbone-induced stability of the system displays itself in form of a temperature-reduced velocity distribution. In contrast to the clustering property, the distribution is determined by the parameter k n.n. in powers of ten, k n.n. = 10 s . Such a property has not been examined systematically, however, for systems which undergo aggregation 48,49 , partly because the latter processes are often carried out under external thermal controls 49 . In this study, we look into the underlying origin of the stability previously found for those systems with s A = −1 and then extend the analysis to systems with s A = −4 over the range of temperatures marked by red dashed line on the left of Fig. 1(e), where the systems undergo processes of aggregation. All simulations are carried out with no controls of temperature.
In conjunction with a macroscopic scenario, the time evolution of a system in molecular dynamics (MD) simulation, consists of a series of step-by-step small changes involve the change δU in internal energy, the macroscopic work δW done by the system and the adsorbed heat δQ at a macroscopic level, in terms of the equality (of the first law) In contrast to the simulation of a fixed equilibrium state, in which the macroscopic work 53,54 and the heat 54,55 are often controlled externally to make those step-by-step changes fluctuating in the vicinity of one single macroscopic state, the series of changes over a non equilibrium relaxation process involves states which may have fundamentally different spatial symmetry or heterogeneity. In MD simulation of equilibrium states, the external controls are based on equilibrium statistical mechanical schemes [53][54][55][56] , designed to generate well defined probability functions 54 , often with controlling designed in accord with the structural features 53,57,58 . In simulating non equilibrium processes with time-varying symmetry or heterogeneity, it is necessary to be cautious in using any external control and better to allow the probability functions be determined by underlying micromechanics, when the interested phenomena, such as the clustering of polymer chains caused by thermodynamic instability, are over limited characteristic length and time scales. The dynamic features dominated by the backbones can be totally preserved if the simulations be carried out without imposing external controls. Though such a strategy has  (c) a snapshot of backbone coordinates from simulation for a system with the angle potential strength K A = 10 −1 ; (d) stages of aggregation: stage I for local lining-ups between segments of two chains, stage II for longitudinal extensions of ordered segments along the chains, stages III and IV for transverse growths to have belt-like two-dimensional local orderings and bundlelike three-dimensional structures, respectively; (e) schematic diagram (in a log-log plot) for equilibrium states (in white color) and steady phases subject to clustering (in grey or yellow color), under the specifically selected setup on the potential U (see Methods). The diagram in plot (e) covers a range of two decades in temperature and a range more than three decades in K A . In simulation, the stages I to IV for clustering are quantitatively identified by the changes in the values of parameters, R 0 , R 1 and R 2 , respectively, over time (see Methods). In Fig. 6, we show the simulation data for systems evolving in time along the red dashed line on the left side of plot (e). In the colored region below the dashed black line, the systems are subject to thermodynamic instability, which drives the initially random and uniform configurations of chains to accumulate. Only those in the yellow region can extend the ordered segments along the backbones (stage II) sufficiently to trigger the formation of transverse domains (stage III) and succeed in aggregation to reach the stage IV with bundled structures in major clusters 48,49 . For the clarity of presentation, the beads shown in the schematic plot (d) do not reflect the sizes of the monomers. The diameter σ (see Methods) of the monomers used in simulations is about three times of the bond length (see Methods). In simulation, the assignment s A = −∞ leads to the cases (not shown in plot (e)) in which we turn off the angle potentials. We also allow to have s = ∞, in which cases the bonds are rigid in fixed lengths. In general, the locations of boundaries between the three regions of plot (e) (in white, grey and yellow colors, respectively) are not sensitive to the strength parameter k n.n. . The snapshot in plot (c) has s = 4 and is at a state marked by 'x' in plot (e). the drawback that the heat δQ contributed by numerical noises is not controlled and the temperature of the system is floating, we do in our studies 51, 52 revealed the stability implicitly supported by the backbone connectivity.
In contrast to the schemes of Monte Carlo simulation 59,60 , where one considers the probability function of microscopic states as a function of total potential energy and temperature T is used as a control parameter, the micromechanics-based schemes of MD simulations allow to measure T by realizing it as a kinetic parameter, determined via the probability density function G(K 1 ) of the thermal kinetic energy K 1 of one particle. The realization can be effectively applied in considering a non equilibrium relaxation process, in which G is a time varying function, specifying the microscopic state of the system at any instant. It is not difficult to realize that the thermodynamic concept of 'quasi-static process' justifies using one single probability function  G for a slow changing non equilibrium process composed of a series of nearly equilibrium states with different values in T, as soon as the kinetic energy being scaled by T.
For an equilibrium state, satisfying Maxwell-Boltzmann statistics, the scaled function  G reads for ν degrees of freedom (k B : the Boltzmann constant) and 1 1 is the Jacobian factor with normalization constant B. We emphasize that the property holds for the thermal part of kinetic energy, on the basis of an isotropic distribution in directions of motions. In a flow, for example, K 1 is contributed by the particle velocity with the local flow velocity being subtracted off.
In this study, we explored the possibility of an extension of such a scheme, to include the situations where backbones dominate the motions of the monomers, in 'softly confined' systems of stiff polymer chains. In these systems, the restricted size and limited time span of the interested phenomena prevent the system from being effectively comprehended by equilibrium statistical mechanics. Instead, from the results of MD simulations for systems with s A = −1 51, 52 , where no external controls are imposed, we have seen some evidences that, the presence of backbone-supported dynamic stability distinguishes the states of reaching 'quasi-steady' situations from those staying transient. We have an extended form for the scaled function Eq. (2), ,1 and the exponential function e x in Eq. (2) replaced by the generalized exponential function e q x , defined by ≡ + − − e q x (1 (1 ) ) q x q 1 1 61,62 . q ranges ≤ ≤ q 1 7 5 and the case q = 1 resumes the ordinary exponential function. In a 'quasi-steady' situation, the values of q for all instants fall in a narrow range around a fixed value q*. It highlights the notion that those states are considered as an integrated situation at non equilibrium.
To be specific, we consider the processes of thermal equilibration between polymer chains and background fluid in six systems with the angle potential by the parameter s A = −1 and with, respectively, the parameter for bonding strengths s = 0, 1, 2, 3, 4 51 and ∞ ref. 52. Over the processes along the red dashed line on the right side  Fig. 1(e), the probability density function P GMB (v) for monomer speed v at each instant is well fitted by the generalized Maxwell-Boltzmann (GMB) speed distribution where the factor j G (v) = 4πv 2 has been used in previous studies 51,52 and A q is the normalization factor and all monomers are assumed to have equal masses, m = 1. Note that, the choice of factor 4πv 2 is compatible with the isotropic Maxwell's distribution of speed (subscript 'IMB' for istropic Maxwell-Botzmann), which is as the special case of q = 1. In this case, the temperature T is determined by . Specifically, we have the number of degrees of freedom ν = 3, for systems of finite s. The factor j G (v) will be discussed further in the last paragraph of Discussion. The "quasi-steady" refers to the situations that the fitted values of q for all instants of a period of time fall over a narrow range measured by ±Δq around a fixed value q*. The distributions for all instants are then well approximated by one single master curve given by refs 51, 52 , 1 . Each of the six systems has N P = 40 chains of n = 100 monomers in length, mixed with a background vapor of N F = 6000 spherical atoms. All of these systems had evolved in time from initial microstates of the same quenched configurations and undergo relaxations toward thermal equilibration between polymer and fluid, which is signaled by the reaching of constancy in the ratio where we put ν P = 3 and ν F = 3, to underscore the numbers of degrees of freedom of the monomer and the fluid atom, respectively. In equilibration, Γ is allowed to reach a value ranging between 201 300 and 1. The lower bound 201 300 is for the cases of rigid bonding (s = ∞), where 99 bond constraints are present for a chain composed of n = 100 monomers. While the evaluation does not take the collective modes as the basis in counting the numbers of degrees of freedom, the constancy in Γ does behave as a convenient signature to identify the reaching of the peculiar equilibration between the two parts, in spite of the varying temperatures for both parts over time. The thermal equilibration within the collection of polymer chains in the mixture, on the other hand, is signaled by the afore-mentioned scaling in the distribution of monomer velocities (Eq. (6)). Specifically, the fitted q's are in a range estimated by q* ± Δq, which are 1.032 ± 0.019 for s = 0, 1.025 ± 0.017 for s = 1, 1.028 ± 0.020 for s = 2, 1.075 ± 0.026 for s = 3, 1.140 ± 0.026 for s = 4 and 1.169 ± 0.047 for s = ∞, respectively. They show an increasing trend with s. Those averages are taken on snapshots over time periods of 30-40τ, in intervals of 0.25τ (see Methods, for the definition of time unit τ).

Figure 2(a) shows how the backbones affect the motion of monomers, in the distribution of the cosine
along the backbone with bond vector → b at that monomer site ( Fig. 1(a)), for those six systems with s A = −1. The enhanced probability emerges in the central part of the distribution when the stiffness increases from s = 0 to s = ∞, indicating the trend in dividing kinetic energy with respect to the backbone coordinates, with a preference in the transverse direction. One important feature of the distribution function ⋅f v b ( ) is the symmetry between ⋅v b and − ⋅v b, as the result of the governance of backbones. It suggests the curved coordinates along the backbone ( Fig. 1(a)) be the proper reference to quantify the motion. Figure 2(b) shows the distribution of the cosine ⋅v b for systems of fixed s A = −4 and s = 2, 3, 4, and ∞, respectively. Again, as s increases, the distribution has larger values near ⋅ =v b 0. One immediate refinement of the information on the anisotropy in motion in reference to the backbone coordinates is to decompose the velocities of individual monomers into the tangential and the transverse components and fit the distributions approximately by one-dimension Maxwell-Boltzmann functions. While there are recognizable deviations between the data and the fitted curves, for those cases with strong anisotropy in motion, such as the systems with s = 4 and s = ∞, in treating them as if they were statistically mutually independent components, the combination into an elliptical Gaussian distribution 63 does provide an acceptable approximation for velocity distribution in three-dimension. Figure 3 shows such fittings for the six systems with s A = −1 by using the 'elliptical Maxwell-Boltzmann' (EMB) distribution, defined by Scientific RepoRts | 7: 3105 | DOI:10.1038/s41598-017-03136-7 where the Jacobian factor and the normalization constant The error function erf(.) is defined by  . While all three fittings provide approximately the same curve as the simulation data for plots (a-c), both the fittings by GMB and by EMB agree with the simulation data better than it is for IMB, for plots (d-f). Note that the values for q shown in the plots are obtained by the fittings over 0.002τ and fall in the ranges q* ± Δq (listed in the text below Eq. (7)) obtained by averaging over time periods of 30-40τ. parallel (marked by subscript '‖') and in two perpendicular directions (marked by subscripts '⊥1' and '⊥2' , respectively) in reference to the curved backbone coordinates, corresponding to b , û and ŵ, respectively ( Fig. 1(a)). In the limit → ⊥ T T , we have j E (v) → 4πv 2 and Eq. (8) converges to the isotropic case, Eq. (5). The values listed for the fitted curves in Fig. 3 show that, for those systems with ⋅f v b ( ) (Fig. 2(a)) deviating from the uniform distribution, their T ⊥1 and T ⊥2 are significantly larger than T ‖ (Fig. 3(d-f)). The polymer chains, under the dominance of backbones, are internally in non-equilibrium situations.
The fitting by EMB (Eq. (8)) is certainly better than that by the isotropic Maxwell-Boltzmann distribution (IMB) (Eq. (5)) and is comparable with the result by Eq. (4). Equation (8) suggests a scenario of thermal motion in a 'tube'-like string 37,[64][65][66] . Note that, while both Eqs (4) and (8) are two-parameter fitting functions, the analysis using Eq. (4) does not require the knowledge of backbone frames and provide a convenient way to reveal the stability of the system, by simply fitting to obtain a scaled distribution function covering all instants of time. Besides, the validity of expression Eq. (4) 67 does not rely on the statistical independence among the components of velocity, along the longitudinal and in the transverse directions. In previous studies 51, 52 , it has been shown that the fitting by Eq. (4) and, thereof, by Eq. (6) for sequences of microscopic states, do identify effectively the non-equilibrium 'quasi-steady' states, generated by untamed numerical noises. In this study, we retained the same  Fig. 1(b)) for l = 1, 2, 4, 8 and 16, observed in the curved coordinates, for mixture systems described in Fig. 2 for the fully isotropic situation. We also list the values for the corresponding q-parameter obtained by the fittings of the generalized Maxwell-Boltzmann (GMB) speed distributions (Eq. (4)).
Our next task is then to extend the temperature-reduced approach to resolve the following issues. Firstly, while the parameters Γ and q provides convenient signatures to find, respectively, the external and internal thermal equilibrations for a collection of polymer chains, we are lacking of knowledge on the underlying origins of the stability leading to the equilibration. Are they from the same origin or not? Secondly, the equilibration reached in six systems considered so far are derived from the same initial preparation. Over longer time scales for sufficient relaxations, allowing to erase the residual memory of initially excited collective modes 68 , in what form will the equipartition of energy prevail in these systems? Thirdly, the non-equilibrium situations considered so far are in forms of temperature varying caused by external source (heat from numerical noises). Whether and how the governance by backbone could sustain the stability under the intrinsic changes in a process of aggregation is subject to be examined. Fourthly, it is not clear whether there is a system in the class of model chains under our consideration does not have a backbone dominant 'quasi-steady' stability?
We proceed to analyze the internal correlations along the chains, to underscore the effects of backbone connectivity. The quantities of which the pair correlation to be analyzed should not be directly dependent on temperature and should reflect the orientations of local coordinates of the curved frame along the backbones.

Correlation-anti-Correlation Symmetry in directions of motions.
We consider the pair correlation between the directions of motion for monomers along a chain, measured by the probability density f l of direction  ( Fig. 4) and q for systems with their parameters (s A , s) marked in the plot, the inner product is evaluated in terms of backbone coordinates, defined in Fig. 1(b); (b) the 2D plot of the mean ⋅ +v v i i 16 evaluated in terms of simulation-box ( Fig. 1(c)) coordinates and q; (c) the normed square deviation  ∞), nevertheless, do show the consistent trend that support the above-mentioned conjecture. The directions of motion in all these eight systems are, therefore, in backbone-constrained isotropy, including the system with (s A , s) = (−4, 3), even though the latter has varying q* over time ( Fig. 6(b)). We also include, in plots (a,b), three data points (red diamonds) for a system with (s A , s) = (−∞, ∞) (i.e. with rigid bond and vanishing angle potentials) and in absence of background fluid 48  tionship between v i and the local orientations at site i ( Fig. 1(a)) and revive v i by resuming the relationship in the new local coordinates at site i + l, as Fig. 1(b)). We briefly denote the inner product  Figure 4 shows the typical curves of f l 's, for l = 1, 2, 4, 8 and 16, computed over time intervals of the same duration as that considered in Fig. 2(a), which is 125 times of the period employed to compute the speed distributions in Fig. 3. The distributions for smaller l's, in all of the six systems, are highly asymmetric between the positive and the negative sides of . indicates the tendency of being parallel in motion for a pair of near neighbors. The tendency is dominantly stronger for the three cases with stronger bonds (Fig. 4(d-f)), in comparing with the other three ( Fig. 4(a-c)), where f l 's are close to be uniform for all l's. The deviation from a uniform distribution increases with the enhanced bonding (by increasing s). For larger l's, the asymmetry fades away, that we see the emergence of symmetrized f l 's, with equal probability to be parallel-bound with a value = > We emphasize that such a symmetry is present only when the correlation is measured in referring to the backbone coordinates ( Fig. 1(a-c)). It displays the isotropy of thermal motions of monomers under the constraint of backbones. Quantitatively, the vanishing mean i i l l 1 1 as l → ∞, is a necessary condition for such a symmetry. For those cases considered in Fig. 4, the magnitudes of their means for a pair separated by l = 16 bonds are less than 0.01 (filled black circles in Fig. 5(a)). The same computation with the inner product defined in the ordinary way for the first moment ⋅ +v v i i l (filled black circles in Fig. 5(b)), by using the simulation-box coordinates, leads to a larger value for each snapshot, than its counterpart In each case, the top panel shows the instantaneous temperatures ⁎ T P and ⁎ T F for polymer and for background fluid, respectively; the middle panel shows the effective number of degrees of freedom ν eff per monomer and the fitted q*; the bottom panel shows the R-parameters: R 0 , R 1 , and R 2 . The colored dashed lines in each bottom panel are plotted to show the reaching of steady values for the R 0 , R 1 and R 2 parameters, respectively. The systems were prepared as uniform configurations, with R 0 , R 1 and R 2 staying close to unity (in t < 0, not shown in the plots). The initiation of the growth is signaled by a rapid falling in R 0 , corresponding to the entry into stage II on the triggering of paralleled segments in stage I ( Fig. 1(d)). It is followed sequentially by the fallings of R 1 and R 2 , which mark the reaching of stage III and stage IV, respectively. The latter (stage IV) is featured by the formation of bundled structures (Fig. 7) which may start immediately after the initiation of the previous one (stage III) 48,49 ; the system in plot (c) has an apparent well separated initiations for the two stages. Each data point for q* is an average on 100 data points over a shift window of Δt = 25τ, which is 100 times of that used for Figs 2(b) and 5, in computing the values of the parameter. q* and ν eff are plotted in spacing of 0.25τ. The curves of both quantities are guided by horizontal dashed lines to show their convergence towards constants, except that of q* in plot (b) for the system with s = 3. In the latter system, q* decreases with time t, in contrast to ν eff which barely changes. Data in an extended simulation of the system are also plotted on the right of in plot (b), which shows continuing decreasing trend in q*. The factor ν eff , on the other hand, encounters a major change separately from that in q*.
Scientific RepoRts | 7: 3105 | DOI:10.1038/s41598-017-03136-7 Fig. 5(a)) in referring to backbone coordinates. The results indicate the unique status of the backbone coordinates in quantifying the correlation properties. The difference between f 16 and the distribution for unconstrained (full) isotropy (abbreviated as 'UI'), f UI (x) = 0.5, measured by the normed square deviation  Fig. 5(c)), indicating the anisotropy in motion is dominated by chain stiffness. Note that, the properties of f l 's basically have no direct dependence on the static properties of the conformation of backbone. The statement is indeed true for the decay length l 0 (the persistent length) of the correlation between the tangential backbone directions ⋅ +b b k k l at the site-indexed distance l, which is an exponentially decaying function in l. While the shapes of distributions f l 's are qualitatively dispersed over a wide range among those cases, the persistent lengths are robustly around l 0 ≈ 10. The systems under quasi-steady situations are, therefore, maintaining a kind of dynamic balance, characterized by the 'parity symmetry' in f between ⋅v b and − ⋅v b (Fig. 2(a)) and by the 'correlation-anti-correlation symmetry' in f l , for large l's, between  Fig. 2(a). It is worth to mention at this point that, in all cases above, the distributions of the tangent directions b of the backbones are isotropic in space. In systems of chains undergoing clustering, the external isotropy become destroyed during the processes. One important issue is then how the validity of Eq. (4) responds to such a change. In the following, we proceed to probe the issues by considering aggregation processes over the whole time spans of relaxations.

Symmetry in Conformation and Isotropy of Motion.
We consider first a case of aggregation process, over which the system is totally out of the status of quasi-steady situation. The process occurs in a pure system of homopolymer chains 48 , with rigid bonding (i.e. s = ∞) and vanishing angle potentials (i.e. s A = −∞), in absence of background fluid (N F = 0). Over the process, the system has neither a well-defined velocity distribution specified by Eq. (4) with a virtually fixed q for monomers; nor has it the symmetry in the distribution f l for direction cosine + v v i i l between a pair of monomers at large separation l along a chain. In Fig. 5(a), the first moment of f l , Eq. (9), for l = 16 (red diamonds), in snapshots of three arbitrarily picked time spots during the process, which possess correspondingly well-separated values in fitted q's, are much larger than those for the other systems plotted in the figure. The result indicates the motions of the monomers are not in 'backbone-constrained isotropy' over the process, during which the backbone tangent orientation b is loosing isotropy.
The observation does not, however, imply that a non-equilibrium system with emerging anisotropy in conformation would necessarily loose the symmetry in The aggregation process occurred in a mixture system of chains of rigidly bonded (s = ∞) monomers, with the tiny but non-vanishing conformation potentials specified by s A = −4 ( Fig. 1(e)), in presence of a background fluid 48 , does display the 'backbone-constrained isotropy' in motion. It possesses virtually fixed q over the process. The data in Fig. 5 (black open circles) display that the first moment of f l for l = 16 is also vanishingly small ( Fig. 5(a)), which indicates the presence of 'backbone-constrained isotropy' in the system.
We extend the analysis to other systems with s A = −4 and consider the cases with s = 2, 3 and 4, respectively. These systems were prepared by the same procedure as those considered in Figs 3 and 4, with polymer chains mixed with background fluid. We monitor the effective number of degrees of freedom per monomer, as the alternative parameter for the ratio Γ (Eq. (7)), in identifying the reaching of thermal equilibration with the background fluid. We would have the convergent value ν eff = ν F = 3 for a collection of non bonded monomers and ν = + eff n n 2 1 for rigid polymer chains, each composed of n monomers. All of the three systems undergo aggregation processes and the mechanisms are quantified by the set of R-parameters (see Methods), where the decreasing in the values of R 0 , R 1 and R 2 48, 49 signals the entering of stages II, III and IV, respectively, depicted in Fig. 1(d). Over the processes (Fig. 6, with typical snapshots shown in Fig. 7), the values of q*, averaged over a shifting time window of width Δt = 25τ, for the systems with s = 2 and 4, respectively, remain constants (middle panels of Fig. 6(a,c) respectively.). The system with s = 3, however, has a decreasing q* (middle panel of Fig. 6(b)), despite of the reaching of equilibration between polymer and background fluid signaled by the constancy in ν eff . An extended simulation (on the right of Fig. 6(b)) shows that the adjustment in the effective degrees of freedom of the polymer chains proceeds as well.
In Figs 2(b) and 5(a-c), we include the data for snapshots at some time spots of these three systems, one snapshot for each of the two systems, with s = 2 and s = 4, respectively (in black open circles) and three snapshots for the system with s = 3 (in blue open circles). Over the aggregation processes, the motions of the monomers along the chains in these three systems stay at the status of having 'backbone-constrained isotropy' , manifested by the presence 'correlation-anti-correlation symmetry' in 16 are all vanishingly small (Fig. 5(a)). While the normed deviation for the snapshots of the systems with s = 2 and 4, respectively, are in agreement with the curve determined by the family with s A = −1 (filled black circles), those for the three snapshots of the system with s = 3 are apparently not on the curve (Fig. 5(c)).
To comprehend the seemingly dispersed situations around the fittings to Eq. (4) among those three systems, which maintain the backbone-constrained isotropy in motions of monomers over the symmetry-breaking non-equilibrium aggregation processes, we go back to examine the implicit assumptions made behind the equations of fittings. Our observations suggest, in the quasi-steady situation under the constraints of backbones, the orientation dependence of velocity is on the angle between v and b , specified by ⋅v b, that is, where P c is a conditional probability for v under the given ⋅v b and ⋅f v b ( ) is the distribution of the variable ⋅v b. The expression Eq. (12) is valid under the condition that the dependence of the joint probability of the three variables v, v and b is independent of the variable b . The condition considers the backbone coordinates as the relevant variables for statistical quantities and is supported by the unique status of the coordinates to demonstrate symmetry between the correlation and the anti-correlation of monomer motions.
Based on Eq. (12) the data analyzed in Fig. 6 suggest the following scenarios. For uniform distributed ⋅v b, ) 0 5, the variable ⋅v b has virtually no dominance on the speed variable v. The distribution P(v) is simply the standard Maxwell's speed distribution for isotropic motions, Eq. (5). Increasing s results in the enhanced deviation from the uniform distribution for ⋅f v b ( ) (Fig. 2). For the two cases with (s A , s) = (−4, 3) and (s A , s) = (−4, 4), respectively, their values of q parameter are greater than unity. The large stiffness in the latter system (s = 4) seems render the conditional probability unaffected by the supposedly enhanced interplays between different chains, in contrast to that for the former system (s = 3), where the q parameter in Eq. (4) varies with the changes with time. Since the decreasing q* is also observed for the system with (s A , s) = (−1, 3), the chains of which do not aggregate, in extending the time span covered in the previous simulation 51 to the same orders considered in Fig. 6(c), the tendency should be mainly due to the reaching of time scales sufficiently to allow the softening of the biased modes that had kept the chains staying in the balanced non-equilibrium situations. Such an internal property for individual chains may less sensitive to the external changes in symmetry of conformation encountered over the aggregation process. The data of the extended evolution in Fig. 6(b) suggest that the time span required for the full relaxation in reaching equilibrium for the system with the set of  Fig. 6 for each of these snapshots. parameters (s A , s) = (−4, 3) would be well beyond the time covered by the current simulation. Such a time of relaxation is even increased for systems with larger s.
To sum up the results of our analysis, the implicitly balance maintained dynamically by backbones can be present for chains with non-vanishing angle potentials (K A > 0 or ≠ −∞ s A ), which underscores the importance of three-body and four-body interactions (described by the bending angle and torsion angle potentials, respectively). For these systems, the chain stiffness determines the statistical properties of the motions of the monomers under the balance. For chains with weak stiffness (smaller s), the thermal motions of monomers follow the standard Maxwell-Boltzmann statistics. In enhancing the chain stiffness (by increasing s), the monomers may stay in biased modes with anisotropic motions, with respect to the backbone coordinates, for prolonged time spans before these biased modes could relax to partition the energy to all allowed modes. We reveal from the simulation data novel forms of symmetry and isotropy for monomer motions under the constraints of backbones for systems either possess or never undergo aggregation processes. Our observations may underlie the empirically observed statistical properties for the speed distribution in non-equilibrium.

Discussion
A thermodynamic quasi-static process involves small changes from one equilibrium state to another, where the displacements at microscopic level can be divided into the contributions producing a macroscopic work δW and those generating uncontrolled heat δQ (Eq. (1)). The latter are distinguished from the former by displaying spatial isotropy in the motions of particles, supported by the uncorrelated randomness. Such a scenario is not necessarily effective for processes occurring in a system of polymer chains with strong backbone-constraints on the motions of the monomers in a spatially restrictive environment with changes occurring over a time span not allowing the system sufficiently probing over entire phase space. The analysis in this study shows, nevertheless, the using of the major language of equilibrium statistical mechanics to describe those non-equilibrium processes is feasible. The scenario of a 'quasi-steady' situation is identified for the system, where the monomer thermal motion characterized by a kind of uncorrelated randomness under the constraints of backbones. Such a property is far beyond triviality in that, the uncorrelated randomness behind the spatial isotropy in thermal motion is identified, in this context, via the correlation-anti-correlation symmetry of directions of motion measured in the backbone coordinates. The symmetry is valid strictly only for directions of motion recorded in the manner referred to the curved backbone coordinates (Fig. 1). The presence of such a property is closely related to the support provided by the three-body and four-body interactions along the backbones, because the system fails to maintain the required implicit stability if the polymer chains have zero angle potentials.
In summary, we have extended the scenario of a 'quasi-steady' non-equilibrium situation 51 , to comprise the implicitly balanced situations over aggregation processes. We have shown in our model studies that backbone connectivity guides not only the symmetry for the conformation of the polymer chains but also the allocation of kinetic energy and, thereof, the thermal statistics. Under the governance of backbones in a softly confined environment, the kinetic energy of the monomers may be trapped in the transverse modes and stay for an extended time span before it able to be partitioned to all modes. Over the time span, a kind of stability is implicitly sustained, leading to the correlation-anti-correlation symmetry in the directions of monomer motions and the emerging of a class of generalized Maxwell-Boltzmann statistics. In the strong stiffness limit, the statistics can be described consistently by one single temperature-reduced scaling distribution (Eq. (3)) over the entire simulation time after initial transients. The same is true for the weak stiffness limit, where the scaling distribution resumes the unit-less Gaussi an form (Eq. (2)). In between the two extremes, we can still have scaled distributions prevailing over shorter time spans. The tendency is independent ofwhether, or not, the polymer chains undergo symmetry-breaking aggregation processes. Figure 3 shows that for systems with strong springs to connect neighboring monomers in polymer chains, the speed distribution of monomers in the parallel (along the direction of the spring to connect two neighboring monomers) and perpendicular directions have different effective temperatures and such systems are in non-equilibrium states. The non-equilibrium behaviors of peptide chains have also been observed in native collagen fibril 69,70 and oxyhemoglobin crystals 71 . When one increase the temperature of native collagen fibril, the measured curves for Young's modulus and logarithmic decrement are very different for different rates to increase the temperature, i.e. the native collagen fibril shows glassy behavior, similar to the spin glass 72 which has slow relaxation; such slow relaxation at low temperatures can be slower than critical slowing down at the critical point 2,15 . Taking Fig. 3 and experimental data of native collagen fibril 69,70 and oxyhemoglobin crystals 71 together, one would like to propose that interacting polymer chains under some situations can help biological systems to maintain in non-equilibrium states 2, 15 , even for very long times such as that of ancient date seeds 13,14 .
In the present study, the total number of monomers and solvent particles is only 10 4 . By using a more powerful computing system and parallel algorithm, Komatsu, et al. 73 had used MD to simulate fluid turbulence with particle number up to 3.779136 × 10 8 . It would be interesting to use computing resources of the order of ref. 73 to simulate polymer-solvent systems with larger numbers of monomers and solvent particles, and with refined structures. One can then study, for example, the influence of polymer chain length and orientation varieties at smaller scales, to underscore their effects on the phase boundary of Fig. 1(e), the speed distribution function of Fig. 3, etc. The effectiveness of the fitting of velocity distributions by Eq. (4) and, thereof, its temperature-reduced form, Eq. (6) suggests the connection between the stability maintained by the backbones revealed in this study and the q-Gaussian based statistical mechanics. Indeed, the q-Gaussian statistics has been found for thermal fluctuations in spatially heterogeneous systems 74,75 , as the generalization of Maxwell-Boltzmann statistics, prevailing in the non-ergodic situations 76 . In some what different context, Brito, et al. used q-exponential function e q x , to study the statistics of degree distribution in classes of spatially constructed complex networks, similar to our observations, the value for the parameter q has been shown to be continually tuned by some strength-related parameter 77 .
Simulations with enhanced data statistics will reveal more insights on issues, such as how the q-Gaussian distribution Eq. (4) of total speed v, is related to the distributions of those for the components v ‖ , v ⊥1 and v ⊥2 . A theoretical construction to produce a q-Gaussian analog of Eq. (8) will clarify more precisely why Eq. (4) is still fine in the cases when the components v ‖ , v ⊥1 and v ⊥2 have different effective temperatures. Note that there are built-in intrinsic correlations among the components in q-Gaussian statistics 67 , in contrast to the scenario behind the elliptical Maxwell-Boltzmann distribution, Eq. (8), where the components are statistically independent and the factor j E (v) is purely geometrical. On account of the algebraic property of q-exponential that x e n q x can be decomposed as the sum of (n + 1) terms each in form of ′ e q ax , the factor j G (v) in Eq. (4) carries effectively a spectrum of q's in the distribution and is not necessary merely a geometric factor. With intensive simulations of larger size systems, one can explore further the statistical properties of correlations by treating the correlations as variables (Figs 2 and 4). It will be helpful to construct the q-generalized statistical relationship between the direction components of velocities and the speed for the anisotropic situations and to refine the knowledge of the statistics behind Fig. 5(c).

Methods
Interaction Potential. The total interaction potential (Fig. 1) of a collection of polymer chains is modeled by In previous work, we separate strength parameters, denoted by s b and s t , respectively, for the bending angle potential and the torsion angle potential. Since we always assign the same value for the two parameters, we use one single parameter s A to unify the notations. The first term is the sum of each Van der Waal potential V r ( ) ij VWij specific to the inter-chain or non-nearest-neighbor-intra-chain or monomer-fluid interaction sites i and j, at a distance r ij . This term is the only one has a site-dependent potential V ij VW, , indexed by ij. In the model 51, 52 , Lennard-Jones potential, 6 , is used. We assign different sets of values for (σ ij ,  ij ) to specify the interaction between monomers, that between fluid atoms and the fluid-monomer interaction, denoted by (σ, ), (σ F , F  ) and (σ PF , PF  ), respectively. In simulation, we assign  PF =  F =  and σ PF = σ and σ σ = In monomer-fluid interaction, however, we drop the attractive part and retain only the r −12 -potential. A tiny amount (5%) of impurity monomers are introduced, to have the same interactions as those of fluid atoms. These monomers are the only sources of heterogeneity for the polymer chains. The second terms of U is the sum of every intra-chain potential between each nearest neighbor pair i and j along each chain, with V n.n. and k n.n. , respectively, defining the function form and the strength of the potential, which are common to all such pairs. V n.n. is realized as a harmonic potential = − . . define the functional forms of the potentials. To cover a range of varieties in the major properties of the backbones, revealing their effects on the clustering, we extend the united-atom model of polyethylene (PE), which has a specified value 67.187 o for the parameter θ 0 of the potential Θ and assign, in our model, the value θ 0 = 112.813 o . We use unit-less strength parameters = K 10 s A A and k n.n. = 10 s to cover ranges of magnitudes in powers of ten. All quantities are expressed in terms of the length parameter σ and the energy parameter  of the Lennard-Jones potential between a pair of non-neighboring monomers along a chain (Fig. 1). For model of polyethylene (K A = 1, k n.n. = 1 and the bond length σ → = . ), the values for σ,  and τ take the values σ = 4.3 Å,  = 57 K − k B (k B : Boltzmann constant) and τ = 2.2 picoseconds, respectively. The time step of integration used in the simulations for all data shown in this paper is δt = 2.5 × 10 −5 τ.. Figure 1(a) shows local coordinates at site i, defined by the direction vector = → | → |  Fig. 1(a)) of site i, with velocity → v i and the local coordinates is revived at site i + l, to obtain the unit vector + v ( ) i i l (in light red color in Fig. 1(b)). The inner product Preparation. With the setup, we can easily control the susceptibility of the system to aggregate, simply by tuning the values of K A . All systems in this study contain N P = 40 chains, each composed of n = 100 monomers, mixed with a background vapor with N F = 6000 background spherical atoms, in a cubic simulation box of size L = 33.39σ on each edge. The persistent length of the chains ranges from 6 to15 times of bond lengths indicating stiff backbones and their end-to-end distance is less than half of L and their gyration of radius is around one-tenth of the fully extended chain length. The linear dimension L of the simulation box is not much larger than these characteristic lengths. Under such a circumstance, it is difficult to justify the assumption made for equilibrium statistical mechanics without controversy. The latter applies to the situation that the size of individual molecules be much smaller than the dimension of the whole system. Such a 'softly confined' allocation in simulation, on the other hand, reflects certain features encountered in confined biological environments. Note that the analysis carried out in Fig. 3 takes data from very short periods because of temperature variations for simulated systems under no thermal controls, in contrast to those in Figs 2(a) and 4, where the quantities are derived from direction vectors. The latter vectors carry no direct dependence on temperature and a time period of 125 times of that in Fig. 3 to take averages for each curve. The intervals for Fig. 2 Fig. 1(d)). Stage II specifies the longitudinal growth to have ribbon-like local structures. The orderings in stages III and IV are in transverse directions to form belt-like (two-dimensional) and bundle-like (three-dimensional) structures, respectively. The initiation of the aggregation process is signaled by a rapid falling in R 0 , corresponding to the triggering of parallelized segments. It is followed by the growths in transverse directions. Since the systems under consideration in this study are dilute in polymer concentration, it is insufficient to form extended planar crystal structures 50 and, instead, the chains grow into bundles 48,49 (Fig. 7). It is these parameters that are used to determine the schematic diagram of steady phases in Fig. 1(e). Tables 1,   2 Table 3. Statistical Quantitties.