Calculation of a fluctuating entropic force by phase space sampling
A polymer chain pinned in space exerts a fluctuating force on the pin point in thermal equilibrium. The average of such fluctuating force is well understood from statistical mechanics as an entropic force, but little is known about the underlying force distribution. Here, we introduce two phase space sampling methods that can produce the equilibrium distribution of instantaneous forces exerted by a terminally pinned polymer. In these methods, both the positions and momenta of mass points representing a freely jointed chain are perturbed in accordance with the spatial constraints and the Boltzmann distribution of total energy. The constraint force for each conformation and momentum is calculated using Lagrangian dynamics. Using terminally pinned chains in space and on a surface, we show that the force distribution is highly asymmetric with both tensile and compressive forces. Most importantly, the mean of the distribution, which is equal to the entropic force, is not the most probable force even for long chains. Our work provides insights into the mechanistic origin of entropic forces, and an efficient computational tool for unbiased sampling of the phase space of a constrained system.
According to the second law of thermodynamics, a system tends to seek a higher entropy state. When the increase in entropy can be gained through a change in a spatial coordinate, an effective force emerges. This force, purely entropic in origin, is a universal phenomenon, not specific to the underlying microscopic Hamiltonian of the system. Examples of entropic force are found in all realms of physics such as the Casimir force Kardar and Golestanian (1999), elastic force, and depletion force Leckband and Israelachvili (2001). The concept of entropic force has even been applied to gravity Jacobson (1995); Verlinde (2011).
Polymers provide an excellent model system to study the entropic force both theoretically and experimentally. They can be attached to a surface, confined in a volume, or constrained in a particular shape. Compared to polymers free in solution, geometrically confined polymers have reduced entropy. This reduction in entropy leads to an entropic force, which may influence the rate of detachment from a surface or the escape time of a polymer from confinement. The confinement-mediated force was measured for a chain in a tube Turner et al. (2002), and calculated using semianalytical means and brownian dynamics simulations Prinsen et al. (2009). A simple heuristic argument estimates its magnitude to be , where is the monomer length.
More involved calculations include computation of the full partition function under confinement and its derivative with respect to displacement Zhao et al. (2011). Alternatively, this force can be calculated from the local concentration of mass points near the surface Guo et al. (2005); Odenheimer et al. (2005). For a chain tethered to a hard wall, the force increases with chain length, but saturates at . For a temperature \SI310\kelvin and a monomer length of \SI3.4\angstrom, this corresponds to \SI12.6\pico\newton. This is consistent with \SI13.3\pico\newton value for the pressure integrated over the entire surface found by Bickel et al.Bickel et al. (2001), where a shorter monomer length of \SI3\angstrom is assumed.
These methods are adequate to calculate the entropic force as the mean thermodynamic force conjugate to displacement. However, entropic force itself is induced by a fluctuating quantity Bartolo et al. (2002). Several examples in biology point to the importance of rare events that deviate from the average behavior Jeppesen et al. (2001); Guo et al. (2011). Likewise, force fluctuations may have a significant impact on chemical and biological processes Koslover and Spakowitz (2012). Therefore, finding the full distribution of instantaneous forces or its higher moments can give insights not easily accessible by statistical mechanics.
In this study, we obtain the equilibrium force distribution due to a flexible polymer pinned in space or to a surface (Fig. 1). The key idea of our approach is to make random moves in the phase space of generalized position coordinates and their conjugate momenta and calculate microscopic forces exerted by the system. This approach is different from conventional Markov chain Monte Carlo methods that typically apply moves only in the position space. Because of the universality of entropic force, we expect our method to be applicable to a wide range of dynamic systems.
ii.1 Entropic force prediction
The entropic force can be found from a thermodynamic calculation
where is the free energy, is the partition function, and is the generalized coordinate of interest. By measuring the change in the partition function as the chain is slightly displaced from the pin point, we can compute the entropic force. As shown in Fig. 1(b), when the tether is lengthened, additional conformations become available. For the full space case, the number of additional conformations is simply proportional to the area available to the first mass point. For the half-space case, we count the fraction of chain conformations compatible with the boundary condition Waters and Kim (2013), and calculate the small increase () upon varying by to numerically approximate Eq. 1.
ii.2 Phase space sampling
Monte Carlo sampling
Beginning from some initial conformation and initial generalized momenta, we can generate a set of conformations for a coarse-grained chain in a heat bath. The Hamiltonian of the chain in terms of generalized positions () and momenta () will be given by
where the first term is a potential energy based on chain deformation, and the second term is the kinetic energy expressed with the mass metric . For a freely jointed chain (FJC), is zero. We use the FJC model in this study because we want to obtain forces purely entropic of origin, but not due to physical interactions. The FJC model effectively coarse-grains an arbitrary polymer into Kuhn-length segments. The mass metric can be found from the Jacobian
depends on in non-Cartesian space. We choose and angles of each link measured in the global reference frame as the generalized coordinates (Fig. 1(c)). These are chosen over internal coordinates, as a change to the orientation of one link in the global frame will only affect two columns of the Jacobian matrix.
At each step, one of the momentum coordinates can be perturbed, and the new kinetic energy computed. The trial step can be accepted or rejected according to the Metropolis criterion on the total energy. The position coordinates can also be perturbed, which will require computing the change in kinetic energy (via a change in ) before the move can be accepted or rejected.
For longer chains, we sample points in the phase space using modal velocity decomposition Jain et al. (). The key idea is that microstates specified by both and can be directly sampled from the probability density function () to build the distribution of any microscopic variable. Our method is different from methods categorized as equilibrium sampling Zuckerman (2011) in that we sample both positions and momenta. For a canonical ensemble, is given by
where . In this form, sampling individual variables and is not straightforward because they are coupled through . Matrix diagonalization can be used to separate variables, but is not computationally efficient. Instead, the Cholesky decomposition can be used to express the symmetric, positive mass metric as
where is a lower triangular matrix. By defining modal velocities as , we can express the probability in a more separable form
In this form, modal velocities () can be sampled independently of each other from a Gaussian distribution. However, positions () become coupled through , which is known as the Fixman correction Fixman (1974). Hence, we adopt a hybrid approach where we first sample a microstate using normally distributed modal velocities and weight the microstate-dependent variable by the Fixman term. Calculating the determinant of the dense metric tensor is computationally intensive, however Fixman demonstrated that this is equivalent to the determinant of a smaller, tri-diagonal metric of the constrained coordinates (for details, see Supplemental Information).
ii.3 Computing constraint force from positions and momenta
We consider a flexible chain terminally pinned in space or to a surface. We can use Lagrangian dynamics to calculate the constraint force as a function of positions and velocities of the mass points. Because distances between all adjacent mass points are constant, there are degrees of freedom. We use three Cartesian coordinates for the first mass point, and angles ( and ) for all other mass points of the chain (Fig. 1 (c)). Choosing generalized coordinates this way simplifies subsequent calculations. The Lagrangian is given by
For the first mass point of the chain to be radially constrained at distance , we introduce a Lagrange multiplier
where . Using in the Euler-Lagrange equation with respect to
The first term on the right hand side vanishes as does not depend on the Cartesian coordinates of the first mass point. Carrying out the differentiation on the left hand side produces an expression
The second derivative of these terms must be solved for using the equation of motion for the angular coordinates.
Solving for gives us
We can express this more concisely as the geodesic equation using the Christoffel connection coefficient, and making use of the symmetry in the first term under exchange of and
Substituting this back in produces the Cartesian components of the constraint force
where is the covariant derivative in the constrained subspace.
iii.1 A single point constraint in space
We first considerd a thermally equilibrated flexible chain with one end pinned to a single point in full three dimensional space. By using the FJC model without a potential term, we can focus solely on entropic contributions to the force. We used the Monte Carlo phase space sampling method to build a canonical ensemble of chains. As shown in Fig. 2(a), kinetic energy is not equally partitioned between the mass points of the chain due to the constraints. The mass point nearest to the pin point has % less kinetic energy whereas the endmost mass point has % more. Most mass points in the middle have on average, which would be expected for two degrees of freedom. Using these chain conformations and momenta, we calculated individual constraint forces from Lagrangian mechanics (Eq. 15). The distribution of constraint forces is highly asymmetric, including both positive and negative values with a sharp cusp. The mean constraint force is a negative pulling force, and does not depend on the chain length. Due to high skewness of the distribution, the mean force () is larger in magnitude than the most probable force (). The standard deviation also has little dependence on the chain length (Fig. 2(d)).
We compared this mean constraint force with the entropic force () predicted by statistical mechanics. To seek more conformations, i.e., to increase entropy, the chain would prefer detachment from the pin point. As a result, the chain exerts a radial entropic force on the pin point. If the distance () between the pin point and the first mass point increases by , the number of conformations will increase in proportion to the spherical surface spanned by the mass point. Therefore, the partition function () is proportional to , and the force at monomer distance is given by Eq. 1
This entropic force is identical to the mean constraint force in magnitude, and is also independent of chain length as shown in Fig. 2(c).
Interestingly, the mean constraint force due to the entire chain is equal to the mean centrifugal force exerted by a single mass point. A radially constrained mass point has only two degrees of freedom, and therefore, according to the equipartition theorem, its mean kinetic energy is equal to . The mean of the centrifugal force () is
Since the energy () of a canonical ensemble is Boltzmann-distributed, the centrifugal force is exponentially distributed.
As shown in Fig. 2(b), this exponential distribution (dashed line) can explain much of the distribution of negative constraint forces. However, positive forces and a significant fraction of large negative forces deviate from this distribution. This result suggests that freely jointed mass points beyond the nearest one can have a collective impact on the pin point.
To understand combined force generation by freely jointed mass points, we examined conformations that produce anomalous force values outside the range attributable to a single monomer. As shown in Fig. 3, instances where the force has the positive sign result from sharp bending and large velocities. For the simple case of a trimer in two dimensions, we can derive a condition for compressive force
This corresponds to the case where the centrifugal force exerted by on towards the pin point is more than twice that exerted by on the pin point in the opposite direction.
iii.2 A single point constraint on a plane
We next considered a flexible chain attached to a surface. This confinement geometry arises in biology problems Petrov et al. (2013); Ludington et al. (2015) and polymer applications Milchev et al. (2010). The phase space sampling method applies the same way except that we only accept chain conformations that do not cross the boundary. We computed constraint forces at different chain lengths. As shown in Fig. 4(a) and (d), the distribution of constraint forces is highly asymmetric, similar to the full space case. Most forces are exerted in the outward direction, but a significant fraction of forces are compressive. The mean constraint force increases slightly as a function of chain length up to ten monomer lengths, but quickly reaches a constant value. The mean constraint force plateaus at about , which is \SI25\percent higher than the full space case.
For longer chains (), phase space sampling becomes computationally expensive. We therefore adopted an equilibrium sampling method to calculate constraint forces. This method requires weighting by the Fixman correction. We computed a raw force distribution (uncorrected, ’’ in Fig. 4), and later applied this correction (’’ in Fig. 4). Applying the Fixman correction lowers the average value of the force, as the largest forces are correlated with the largest extensions, which have a low weighting factor. This equilibrium method after correction produces similar distributions to the Monte Carlo method, highly asymmetric with more frequent negative than positive forces.
To compute the corresponding entropic force, we consider the increase in the accepted chain conformations upon increasing the distance to the first mass point. The partition function can then be numerically differentiated using Eq. 16 to obtain the entropic force. The result did not depend on the size of displacement chosen. As shown by the red dashed line in Fig. 4(b) and (e), the entropic forces based on conformational space only are slightly larger than the mean constraint forces ( vs. ).
We computed constraint forces exerted on a FJC in full space or half space by classical and statistical mechanical methods. Statistical mechanics was used to describe thermal fluctuations of the chain microstate, and classical mechanics was used to calculate the force required to constrain the chain against these thermal fluctuations. We found that the mean of constraint forces is equal to the entropic force calculated from the spatial derivative of the free energy. However, because constraint forces are asymmetrically distributed, the mean force is not equal to the most probable force. The force distribution is largely independent of the chain length beyond ten monomer lengths.
In this study, we treat the FJC as a canonical ensemble, which exchanges heat with the environment at a constant temperature. Our equilibration of microstates in the phase space is similar to the Andersen thermostat Andersen (1980) used in molecular dynamics with the added capability of simultaneously satisfying geometric constraints. At any instant, the mass points of the chain move with certain velocities at certain positions under the constraints set by the FJC model and the boundary. At a different time point, the mass points will adopt a different conformation with a different set of velocities as a result from collisions with the surrounding molecules of the heat bath. For each set of positions and momenta reset as a result of heat exchange, we calculate the microscopic constraint force.
The most nontrivial result from this study is that the force distribution is not Gaussian even in the long chain limit. Because the mass points are constrained in the freely jointed chain model, their impact on the tether are neither identical nor independent. As a result, the constraint force does not obey the central limit theorem. Consequently, the most probable value and the average value of force do not coincide with each other. Therefore, the entropic force that we derive from statistical mechanics might not be the most relevant quantity in all circumstances.
The fluctuation or second cumulant of an extensive thermodynamic variable can be easily calculated from the partition function when its conjugate intensive variable is fixed Marconi et al. (2008). Textbook examples are energy fluctuation in a canonical ensemble and particle number fluctuation in a grand canonical ensemble. On the other hand, the fluctuation of an intensive variable such as force cannot be derived from the partition function when the conjugate extensive variable is fixed, i.e. in a fixed-distance ensemble Rudoi and Sukhanov (2000); Gerland et al. (2001). Here, we use the FJC terminally constrained in full space to illustrate this point.
The constrained partition function () is given by (See Supplemental Information.)
where is the Hamiltonian of the constrained chain. The mass metric is a matrix that depends on global angles . The dependence of the partition function on the distance between the end and the pin point () is explicity shown. Entropic force () arises because the constrained partition function increases with , which is reflected by the thermodynamic relation:
The quantity in the constrained ensemble average can be identified with a microscopic force, which we denote as . The second moment of can be shown to be
Whereas the first term can be calculated from the constrained partition function, the second term depends on the microstate of the system. If the rigid constraint is relaxed to a spring, the second term is simply the stiffness () of the spring. In this case, force fluctuations can be easily calculated, similar to those for a particle trapped in a harmonic potential Manosas and Ritort (2005); Mehraeen and Spakowitz (2012). However, the force fluctuation diverges in the stiff limit () Gerland et al. (2003). Therefore, to obtain force fluctuations with a hard constraint, one has to resort to computational means.
We showed that for a terminally pinned chain in full space, our phase space sampling method produces constraint forces whose mean is equal to the entropic force . Using Eq. 20b, the entropic force can be worked out analytically. Integrating out the momentum dependence of , one obtains
depends on through the mass metric . It can be shown that (Supplemental Information), and therefore the entropic force is given by
In comparison, the entropic force due to the pinned chain in a half space is difficult to derive because the limits of integration imposed by the half plane depend on in some complex fashion.
where is the height of the pin point from the surface. By taking the derivative of with respect to , a vertical entropic force can be derived. This vertical force increases with chain length for short chains, but is nearly constant at for long chains. This force, however, is expected to be different from the mean constraint force obtained from our phase sampling methods for several reasons. First, the method of images yields only the vertical component of the entropic force out of mathematical convenience. Second, the analytical expression (Eq. 25) is correct only in the Gaussian limit of long chains. Third, the conformational partition function in Eq. 25 does not include the kinetic effect. However, it is unclear if removing the kinetic effect would significantly alter the mean force prediction. Therefore, we numerically computed the conformational partition function as a function of instead of , and obtained a thermodynamic force conjugate to using Eq. 16. This force is only slightly larger than the entropic force obtained by phase space sampling as shown in Fig. 4(b) and (e).
As the chain length becomes longer, the Monte Carlo method becomes computationally expensive. A faster method is to sample conformations directly from an a priori distribution. For example, to simulate a canonical ensemble of ideal gas, one can randomly sample velocities from the Maxwell-Boltzmann distribution and assign them to randomly positioned particles. Such equilibrium sampling, however, is not straightforward for a constrained system due to the coupling of through . By decomposing , we can obtain modal velocities that individually obey a Gaussian distribution. Such a transformation skews the distribution by the Fixman term, which must be corrected for to unbiasedly sample the phase space. Our results show that this correction is warranted (Fig. 4(e) and (f)).
Both phase space sampling methods we introduce in this work are based on theoretically correct interpretations of the constrained partition function. The Monte Carlo phase space sampling method has the advantage of not requiring the Fixman correction to produce the correct average force prediciton. Additionally, this method may prove more useful for geometries such as polymer chains with a constant end-to-end distance. In such cases, selecting conformations that fit a given distribution is difficult but perturbing a conformation that already lies within that distribution is comparatively easy. The chief disadvantage of this method is the computational time required to perturb the chain, and the large number of steps required to adequately explore the conformational space.
In contrast, the equilibrium phase space sampling, where uncorrelated conformations are selected, and momenta assigned after, is much quicker for free chains. This lack of correlation may more easily lend this method to large-scale parallelization as well. One minor disadvantage of this method is that it introduces a discrepancy in the mean force, unless the sampled conformations are weighted by the Fixman correction.
We introduced two different computational methods to sample microstates of a freely jointed chain in an unbiased manner. Unlike mainstream Monte Carlo methods that explore the conformational space only, our methods explore the complete phase space to obtain full distribution of any velocity-dependent microscopic variable in thermal equilibrium. We applied these methods to a terminally pinned chain in full space and half space to calculate constraint forces. The force distribution was non-Gaussian and asymmetric with both tensile and compressive forces. Most notably, the most likely force was smaller in magnitude than the mean force, which is more commonly considered both theoretically and experimentally. The constraint force exhibited little chain-length dependence beyond monomer lengths. The presence of a plane leads to a larger mean constraint force. This work can be extended to investigate fluctuating forces using more realistic polymer models and confinement.
The authors acknowledge financial support from Georgia Institute of Technology and the Burroughs Wellcome Fund Career Award at the Scientific Interface. We thank Tung Le, Jiyoun Jeong, and the rest of the Kim lab for careful reading of the manuscript. We also thank Dr. Kurt Wiesenfeld and Dr. Toan Nguyen for helpful discussions.
.1 Dependence of the mass metric on the end constraint
The Cartesian coordinates of the n-th mass point are given by ().
It is assumed that all adjacent mass points are joined by an inextensible linker of unit length. The spherical coordinates of the first mass point are given by
We define two sets of generalized coordinates and a subset of unconstrained coordinates where . The Jacobian matrix with respect to can be constructed from smaller Jacobians for individual mass points. The Jacobian matrix is a matrix for the first mass point:
and a matrix for the rest:
The total Jacobian can be expressed as a block matrix:
The mass matric becomes a square matrix:
where -dependence is completely isolated in the first block row and column. Expanding the block matrices, one can obtain -dependence of the matrix elements:
where denotes terms proportional to . We can obtain a submatrix by excluding the first row and column of :
It is straightforward to show , which gives rise to an entropic force upon constraining . To demonstrate this, we transform the partition function for the unconstrained system to the constrained one den Otter (2013). We use the Hamiltonian of a freely jointed chain without the end constraint to express the partition function:
is a vector of momenta conjugate to . Delta functions are introduced to apply the end constraint on :
where -dependence of the constrained partition function is explicitly shown. The delta functions can be handled by decomposing the mass metric and using the transformation rule
Integrating with , we can obtain the -dependence of the partition function:
.2 Calculation of the metric determinant
Chain conformations acquired from equilibirum sampling of the conformational space, without consideration of momentum space, will require a weighting factor given by the determinant of the metric tensor in the constrained space. The metric in the constrained space is generally dense, and so calculating the determinant is computationally intensive. Fixman Fixman (1974) made use of the simple form of the metric determinant in the unconstrained space, to find the desired from the determinant in the smaller, orthogonal subspace spanned by the constrained dimensions.
This is further simplified by considering the inverse of the metric for the orthogonal subspace of constrained dimensions. This has the form
where is the bond length between points and . The terms in the sum are only non-zero when and . This ultimately produces a tridiagonal matrix
Expanding this matrix in minors gives a recursive formula for the determinant , where is the unit relative position vector pointing from (n-1)-th mass point to n-th mass point.
For equal mass points, we begin with the base cases . Note this is slightly different than the Fixman case () as we are also constraining the position of the first mass point. Once we have computed , the determinant is easily found from
- preprint: APS/123-QED
- M. Kardar and R. Golestanian, Reviews of Modern Physics 71, 1233 (1999).
- D. Leckband and J. Israelachvili, Quarterly reviews of biophysics 34, 105 (2001).
- T. Jacobson, Physical Review Letters 75, 1260 (1995).
- E. Verlinde, Journal of High Energy Physics 2011, 1 (2011).
- S. Turner, M. Cabodi, and H. Craighead, Physical Review Letters 88, 128103 (2002).
- P. Prinsen, L. T. Fang, A. M. Yoffe, C. M. Knobler, and W. M. Gelbart, The Journal of Physical Chemistry B 113, 3873 (2009).
- S.-L. Zhao, J. Wu, D. Gao, and J. Wu, The Journal of chemical physics 134, 065103 (2011).
- K. Guo, F. Qiu, H. Zhang, and Y. Yang, The Journal of chemical physics 123, 074906 (2005).
- J. Odenheimer, D. W. Heermann, and M. Brill, International Journal of Modern Physics C 16, 1561 (2005).
- T. Bickel, C. Jeppesen, and C. Marques, The European Physical Journal E 4, 33 (2001).
- D. Bartolo, A. Ajdari, J.-B. Fournier, and R. Golestanian, Physical review letters 89, 230601 (2002).
- C. Jeppesen, J. Y. Wong, T. L. Kuhl, J. N. Israelachvili, N. Mullah, S. Zalipsky, and C. M. Marques, Science 293, 465 (2001).
- Z. Guo, M. Gibson, S. Sitha, S. Chu, and U. Mohanty, Proceedings of the National Academy of Sciences 108, 3947 (2011).
- E. F. Koslover and A. J. Spakowitz, Physical Review E 86, 011906 (2012).
- J. T. Waters and H. D. Kim, Macromolecules 46, 6659 (2013).
- A. Jain, I.-H. Park, and N. Vaidehi, Journal of chemical theory and computation 8, 2581.
- D. M. Zuckerman, Annual Review of Biophysics 40, 41 (2011).
- M. Fixman, Proceedings of the National Academy of Sciences 71, 3050 (1974).
- A. S. Petrov, S. S. Douglas, and S. C. Harvey, Journal of Physics: Condensed Matter 25, 115101 (2013).
- W. B. Ludington, H. Ishikawa, Y. V. Serebrenik, A. Ritter, R. A. Hernandez-Lopez, J. Gunzenhauser, E. Kannegaard, and W. F. Marshall, Biophysical journal 108, 1361 (2015).
- A. Milchev, L. Klushin, A. Skvortsov, and K. Binder, Macromolecules 43, 6877 (2010).
- H. C. Andersen, The Journal of chemical physics 72, 2384 (1980).
- U. M. B. Marconi, A. Puglisi, L. Rondoni, and A. Vulpiani, Physics reports 461, 111 (2008).
- Y. G. Rudoi and A. D. Sukhanov, Physics-Uspekhi 43, 1169 (2000).
- U. Gerland, R. Bundschuh, and T. Hwa, Biophysical Journal 81, 1324 (2001).
- M. Manosas and F. Ritort, Biophysical journal 88, 3224 (2005).
- S. Mehraeen and A. J. Spakowitz, Physical Review E 86, 021902 (2012).
- U. Gerland, R. Bundschuh, and T. Hwa, Biophysical journal 84, 2831 (2003).
- Y. Hammer and Y. Kantor, Physical Review E 89, 022601 (2014).
- W. K. den Otter, Journal of chemical theory and computation 9, 3861 (2013).