Simple flexible polymers in a spherical cage

We report the results of Monte Carlo simulations investigating the effect of a spherical confinement within a simple model for a flexible homopolymer. We use the parallel tempering method combined with multi-histogram reweighting analysis and multicanonical simulations to investigate thermodynamical observables over a broad range of temperatures, which enables us to describe the behavior of the polymer and to locate the freezing and collapse transitions. We find a strong effect of the spherical confinement on the location of the collapse transition, whereas the freezing transition is hardly effected.


Introduction
The behavior of proteins in confinements has been studied in theory and experiments for a while. It is a rewarding topic for research because crowded environments such as caverns, cells, micelles, etc., are the natural habitat of biopolymers, and the structural conformation has quite an impact on important subjects such as building biosensors [1], packaging of DNA [2] or the folding behavior of proteins [3][4][5][6]. In this work we are looking at the behavior of a polymer captured in a steric sphere, which can be considered as a simple model for a polymer in a micelle, chaperonin-like cage or small pore in a synthetic matrix, without a complex thermodynamic behavior of the confining structure itself. There have been simulations with Gõ-like protein models such as β-barrel or β-hairpin proteins and some others [3][4][5][6][7] in similar confinements. To get a general overview on the effects of the confinement, we discard the complexity of 20 different amino acids, which leads to a large variety of realizations for proteins, or an enormous amount of different building blocks for synthetic polymers. Instead, we use a simple beadstick homopolymer model, which gives a good overview on general characteristics. As a first approach, we model the sphere as a steric wall without any attractive or repulsive potential. We monitor the change of the collapse and freezings transitions and their temperatures T Θ c and T F c induced by the reduction of the translational entropy and the available space due to the sphere compared with the free polymer. Although it is a relatively simple model, the energy landscape is complex enough and the density of states ranges over many orders of magnitude. Thus, advanced Monte Carlo techniques are necessary to systematically investigate the thermodynamic behavior of energetic and conformational observables.
The rest of the paper is organized as follows. In section 2 we describe the used model and observables in detail, and in section 3 we briefly review the parallel tempering and multicanonical simulation methods. Afterwards in section 4 we present our results, and in the last section 5 we give a short conclusion. angle between two bonds θ i is defined as cos θ i = ( r i+1 − r i ) · ( r i+2 − r i+1 ). As for lattice models, we neglect any bond vibrations, and adjacent monomers are connected via fixed bonds, the distances between these monomers | r i − r i+1 | being set to unity. The excluded volume and attractive parts of the monomermonomer interaction are modeled by a 12-6 Lennard-Jones potential for all non-adjacent monomers, the stiffness is introduced via a bending potential and the confinement by suppressing any state where at least one monomer is located outside a sphere centered around the origin. Summarizing, the Hamiltonian consists of three terms, with the Lennard-Jones part being of the common 12-6 form where r i j = | r i − r j | is the distance between two monomers. The bending energy is given by the usual cosine potential where the parameter κ enables us to adjust the stiffness of the polymer. In the following simulations we set κ to 0.25, thus the polymer is very flexible. The sphere is modeled via where R S is the radius of the sphere which ranges in our simulation from 2 to 12. The polymer is allowed to move freely inside this sphere. Too small spheres lead to a completely unphysical behavior, because the polymer is pressed into conformations smaller than the crystal conformations, which are reported for a similar model in [11,12], and thus the excluded volume part leads to extremely high energies. This gives a lower bound for R S , the upper limit is chosen so that the behavior of the polymer hardly differs from the bulk behavior. In order to observe the freezing and collapse transition and to describe the conformational behavior depending on the radius of the sphere, we choose the following observables. For the freezing transition, the energetic observables are ideal. We measure both parts of the energy E LJ and E bend separately and, of course, the fluctuations of these quantities, C v = d dT 〈E 〉. Additionally, the squared radius of gyration give a good description of the conformational behavior. The maximum of the heat capacity C v is a good indicator of the freezing transition, because at that temperature the polymer moves into a crystal-like structure, which is associated with a strong energy drop induced by the Lennard-Jones potential. The maxima of d dT 〈R 2 gyr 〉 and d dT 〈R 2 ee 〉 are good indicators for the collapse transition, at which the polymer changes its conformation from an extended form to a globular one.

Simulation methods
Although we consider a simple polymer model, its phase space is so complex that the Metropolis Monte Carlo method will lead to misleading results at low temperatures or near pseudo phase transitions. We use two advanced Monte Carlo methods to cope with this problem. A recent overview of these problems is given in [13]. The first method is parallel tempering Monte Carlo sampling, the principle idea is originally described in [14,15] and the algorithm itself in [16,17]. The second method is multicanonical Monte Carlo (MUCA) sampling [18,19]. We use the first one to get a good overview of a broad temperature range, and MUCA to check our results especially near the first-order like freezing transition and at low temperatures. We will briefly introduce both methods here.

Parallel tempering method
The clue of the parallel tempering method is quite simple. One runs several separate Metropolis Monte Carlo simulations, each of them at a different temperature in parallel. Every now and then two replicas are allowed to exchange their conformation with probability p swap = min 1, e ∆β∆E , (3.1) where ∆β is the difference in the inverse temperature of the replicas and ∆E is the difference in the energy of the replicas. In an implementation of the parallel tempering method, one will not exchange the complete state of the system. Instead, one would exchange just the temperature and do a little bit bookkeeping to get everything done right. This procedure is summarized in figure 1. At the end, the whole simulation yields separate time series for each temperature, which is a good starting point for the multi-histogram reweighting technique (WHAM) [20,21]. From the individual energy histograms at every temperature, WHAM calculates the density of states Ω(E ) in an iterative way. As a good starting point for this iteration, we use Ω(E ) obtained by a direct histogram reweighting method [22] which gives a good first estimate for Ω(E ) and, therefore, leads to a faster convergence. The parallel tempering method benefits from the possibility that a single Metropolis Monte Carlo simulation, which is possibly stuck in a conformation at low temperature or near a phase transition, can exchange its state with a replica from a higher temperature and thereby overcome its stuck state. This exchange is only possible if the energy distributions of the two temperatures have a sufficient overlap, which means that the temperature difference between two neighboring replicas should be small enough. This is also the reason for the weakness of the parallel tempering method at low temperatures and at first-order phase transitions. At low temperatures, the energy histograms become very narrow and one needs many different replicas to cover a broad temperature range. At first-order phase transitions, the energy distribution is double peaked with an extremely suppressed regime between the peaks, so that the parallel tempering method still suffers here from the weakness of the Metropolis Monte Carlo method which has the problem of overcoming this extremely suppressed region, but MUCA is capable of countering exactly this problem.

Multicanonical Monte Carlo sampling
The multicanonical method allows one to use arbitrary configuration weights instead of Boltzmann weights to sample the phase space of the system. Therefore, the canonical partition function is modified: With this modification one can try to adjust the weights W (E ) in such a way that the simulation spends equal amounts of time at each energy. To obtain this, the configuration weight should be equal to the inverse density of states: W (E ) = Ω −1 (E ). The density of states is naturally unknown before the simulation.

43008-3
Therefore, the weights should be somehow calculated during the simulation. A possible way to do this is by iteration. The simplest approach is to start at arbitrary weights, run a simulation with this weight, calculate the energy distribution H (E ) and modify the weights via W (n+1) (E ) = W (n) (E )/H (n) (E ). This procedure is repeated until the resulting histograms are flat and span the desired energy range. At the end, one can reweight the result from an equilibrium production run with the last weights to every temperature whose Boltzmann energy distribution lies within the flat energy histogram. A possible method for this is time-series reweighting, where every measured observable is weighted by W −1 (E )e −βE which results in the following formula:

Results
To locate the pseudo phase transition, we consider the maxima of the temperature derivative of E , R 2 gyr and R 2 ee . For short chain lengths, one can see both pseudo phase transitions in the temperaturederivative of 〈R 2 gyr 〉. For longer chains, the peaks for the freezing transition are suppressed, but still visible in the heat capacity C v , see figure 2. The qualitative behavior of d dT 〈R 2 ee 〉 is the same as that of d dT 〈R 2 gyr 〉.
The effects of the sphere on the elongation of the polymer are easily predictable. In the extended phase, above the collapse transition, the extension of the polymer is clearly limited by the sphere. This effect is still visible but reduced in the collapsed phase, between the freezing and collapse transition, and hardly visible in the frozen phase, see figure 3.
In figure 2, we see that T F,N max , which denotes the temperature of the maximum of the peak at the freezing transition for a fixed N , and the width of the peaks remain similar for different R S , except for very small R S where the polymer is pressed into very narrow states. For these small spheres, the collapse transition vanishes completely, neither the peak of d dT 〈R 2 gyr 〉 nor the shoulder in C v exists. It depends on the length of the polymer at which R S this effect takes place.
On the other hand, the collapse transition and its temperature T Θ,N max are strongly effected by the confinement. The peaks of d dT 〈R 2 gyr 〉 decrease, become broader and shift to lower temperatures as R S decreases. A decreasing radius of the sphere pushes the polymer into more collapsed conformations even above the collapse transition. Thus, the difference in conformational observables between the collapsed and the extended phase decreases, which explains the broader and lower peaks. The direction of the shift of T Θ,N max is opposite to the behavior of models for proteins reported in [3][4][5][6][7], where the folding temperature increases with a decrease of the available space. We have simulated a flexible polymer, whereas these works handle relative short proteins which are more in the semi-flexible or stiff regime. The different stiffness is a possible explanation for the different behavior. To quantitatively study the shift of T Θ,N  (15) . 3.25 , where R 0 is the size of the polymer and L the length of a confining cylinder, has been reported in [5] for certain protein models. In our case, R 0 ∝ N ν holds with ν = 1/2 as for a random walk: At the collapse transition, the polymer becomes extended and thus does not "feel" the self-avoidance, and in three dimensions it effectively acts as a random walker where ν = 1/2 (up to logarithmic corrections).

43008-4
Simple flexible polymers in a spherical cage

Conclusion
We have presented a Monte Carlo study of the effects on the pseudo phase transitions of a flexible polymer caused by a steric confinement. Advanced Monte Carlo techniques are used to get a detailed estimate for the heat capacity and radius of gyration and its temperature derivative. It is found that the confinement has hardly any effect on the behavior and the location of the freezing transition, whereas for the collapse transition, the behavior and the location change significantly. Due to the loss of translational entropy and the reduction of possible extended states, the transition becomes less and less pronounced with a decreasing radius of the confining sphere. We found a scaling law for the shift of the location of the collapse transition, which holds for all simulated polymer lengths. This shift is directed towards lower temperatures with decreasing radius of the sphere, which is opposite to what has been claimed in other works simulating more realistic models of proteins. One possible reason is that these proteins are much stiffer than the polymer we simulated in this work. An intended future task is to find out where this difference comes exactly from.