Full Statistics of Conjugated Thermodynamic Ensembles in Chains of Bistable Units

The statistical mechanics and the thermodynamics of small systems are characterized by the non-equivalence of the statistical ensembles. When concerning a polymer chain or an arbitrary chain of independent units, this concept leads to different force-extension responses for the isotensional (Gibbs) and the isometric (Helmholtz) thermodynamic ensembles for a limited number of units (far from the thermodynamic limit). While the average force-extension response has been largely investigated in both Gibbs and Helmholtz ensembles, the full statistical characterization of this thermo-mechanical behavior has not been approached by evaluating the corresponding probability densities. Therefore, we elaborate in this paper a technique for obtaining the probability density of the extension when force is applied (Gibbs ensemble) and the probability density of the force when the extension is prescribed (Helmholtz ensemble). This methodology, here developed at thermodynamic equilibrium, is applied to a specific chain composed of units characterized by a bistable potential energy, which is able to mimic the folding and unfolding of several macromolecules of biological origin.


Introduction
The recent developments of thermodynamics and statistical mechanics concern the thermodynamics of small systems, kept far from the thermodynamic limit, and the stochastic thermodynamics, which is based on Langevin or stochastic differential equations.In the first theoretical approach, the small size of the system is carefully taken into account in order to analyze its effects on the overall behavior of the system [1,2] and, in particular, on the force-extension response in the case of macromolecular chains.One interesting feature of small systems thermodynamics is the non-equivalence of the ensembles for finite sizes of the system, and the convergence to the equivalence of the ensembles in the thermodynamic limit [3][4][5][6].In the second theoretical approach, the out-of-equilibrium statistical mechanics is introduced by means of the Langevin and Fokker-Planck equations, which represent the stochastic evolution of the phase-space variables and of their probability density, respectively [7][8][9][10].In this context, the first and the second principles of the thermodynamics can be re-demonstrated [11][12][13][14] and other important fluctuation-dissipation theorems have been elaborated [15][16][17][18][19][20][21].These results follow from the pioneering Sekimoto idea of the microscopic heat rate along a Brownian system trajectory [22,23].Concerning the Brownian trajectory of a particle an interesting investigation concerns the generalization of the principle of the least action in a probabilistic situation, which is equivalent to the principle of maximization of uncertainty associated with the stochastic motion [24].Also, the well-known quantum uncertainty relation can be proven to hold for non-quantum but stochastic trajectories of a Brownian particle [25].Furthermore, the entropy generation during the stochastic evolution of a system has been studied by means of the Gouy-Stodola theorem [26] and applied to the biological context to model molecular machines in the orginal way [27] and to study the control and regulation of temperature in cells [28].Other approaches for investigating the behavior of molecular motors are based on the over-damped Langevin equation and have been successfully compared to the experimental data of the F 1 -ATPase motor [29,30].
Typically, the units or elements of polymers and of other macromolecules may exhibit a bistable behavior or not, depending on their internal chemical structure.The general behavior, and in particular, the force-extension response of chains without bistability is nowadays reasonably well understood [4,[51][52][53][54].On the other hand, the complexity of chains with bistable units has been recently revealed through force-spectroscopy experiments and is the subject of promising research efforts.Indeed, the conformational transition between two states of each chain unit has been observed in polypeptides, nucleic acids and other molecules.The possibility of measuring the dynamic response of bistable systems is very important for investigating the out-of-equilibrium statistical mechanics since the coupling and/or the competition between the purely mechanical characteristic times and the chemical characteristic times induced by the barrier separating the two states [55] can be directly probed and compared with theoretical results (see, e.g., [18,56,57]).Interestingly enough, for relatively short bistable molecular chains, the applied boundary conditions play an important role in defining their overall response [58][59][60][61].
The first boundary condition corresponds to experiments conducted at constant applied force.It means that, in this case, one uses soft devices (low values of the intrinsic elastic constant of the devices) and the experiments are called isotensional.This configuration corresponds to the Gibbs statistical ensemble of the statistical mechanics, and leads to a plateau-like force-extension curve.The threshold force related to this plateau must be interpreted by the synchronized transition of all the chain units [46,[62][63][64][65][66][67].The second boundary condition corresponds to experiments conducted at prescribed displacement.This situation can be obtained by hard devices (high values of the intrinsic elastic constant of the devices) and experiments are called isometric.The process represents a realization of the Helmholtz statistical ensemble of the statistical mechanics, and the corresponding force-extension curve shows a sawtooth-like shape.This behavior proves that the units unfold sequentially in response to the increasing extension [39,41,[66][67][68][69][70][71][72][73].In any case, the differences between isotensional and isometric force-extension curves, or equivalently between Gibbs ans Helmholtz ensembles, disappear if the number of units is very large since, in the thermodynamic limit, the Gibbs and Helmholtz ensembles are statistically equivalent, as largely discussed in the recent literature [3,6].
Typically in the theoretical analyses conducted to study the behavior of two-state systems under isotensional or isometric conditions (see previous works), the considered quantities correspond to the average values of the fluctuating variables.It means that one considers the average extension in the Gibbs ensemble and the average force in the Helmholtz ensemble.However, it is important to study the actual distributions of these fluctuating or stochastic variables in order to better understand the random behavior of these systems and to draw more refined comparisons with experiments.Indeed, it is important to underline that the experimental activities above outlined may probe not only the average values of the relevant quantities but also their actual distribution.Basically, this is achieved by a very large statistics of the trajectories of the system under investigation, which allows for a good exploration of the phase space and, consequently, for the determination of the pertinent probability densities.Therefore, in this paper we propose a methodology to determine the exact distributions or probability densities of the pertinent quantities defined in both the Gibbs and the Helmholtz ensembles.This methodology is developed here for systems at thermodynamic equilibrium, as discussed below.In particular, for the Gibbs ensemble, we determine the distribution of the couple ( ẋN , x N ) where x N is the extension of a chain of N bistable elements (under applied deterministic force), and, for the Helmholtz ensemble, the distribution of ( ḟ , f ) where f is the measured force (under prescribed deterministic extension).When the number of the units approaches infinity (thermodynamic limit), the two ensembles become equivalent as previously stated [3,6].This means that the force-extension responses converge to the same curves.Conversely, the above defined probability densities are not the same for N → ∞ since they are defined through different variables and can not be directly compared.The applied method is based on the spin variables approach, recently introduced to deal with bistable or multistable systems at thermodynamic equilibrium.The idea based on the spin variables approach consists in considering a discrete variable (or spin) associated to each bistable unit, able to define the state of the unit itself.It means that an arbitrary bistable potential energy function describing a chain unit can be reasonably approximated by two spring-like quadratic potentials representing the two wells and, therefore, the switching between them is governed by the behavior of an ad-hoc discrete or spin variable [74].Of course, when we adopt the approximation of the energy wells with two quadratic functions, we lose the information about the energy barrier between the wells and therefore we can not use this version of our model to deal with out-of-equilibrium regimes [55].This approach has been recently used to investigate the properties of several two-state systems and macromolecular chains [74][75][76][77][78].Both the Gibbs and the Helmholtz ensembles can be studied by the spin variables methodology, permitting to draw direct comparisons between isotensional and isometric conditions, provided that we work at thermodynamic equilibrium.While the application of this technique to the Gibbs ensemble is more direct since the integration of the partition functions can be typically performed in closed form without particular difficulties, the approach used for the Helmholtz ensemble is more involved.Indeed, in this case the partition function can not be directly integrated but it can be obtained as the Fourier transform of the Gibbs partition function, analytically continued on the complex plane.The mathematical details about this idea can be found in Refs.[74][75][76].In the present analysis, this approach leads to closed form expressions for the probability densities defined above, and the final results can be interpreted by introducing a form of duality between the two ensembles, useful to better understand the specific features of the isotensional and isometric conditions.The system considered in this work is quite simple and it should be viewed as a toy-model best used to introduce the concepts and discuss the results.Of course, this model and its analysis could be generalized by taking into account more refined features (energy difference between the states, heterogeneity, two-or three-dimensionality, cooperativity among the units and out-of-equilibrium evolution) in order to represent more realistic systems.Here, we reduced the complexity as far as possible with the aim of presenting the adopted methodology as effectively as possible.
The structure of the paper is the following.In Section 2, we introduce the investigated model and we determine the partition functions under both the Gibbs and the Helmholtz ensembles.Here, we also discuss the force-extension relation for the two ensembles.In Section 3, we introduce the complete probability density for the system in the whole phase space.This is a preliminary information exploited afterwards to deal with the specific distributions of the two ensembles.Indeed, in Section 4, we obtain the probability density of the couple ( ẋN , x N ) versus f within the Gibbs ensemble, and in Section 5, we get the probability density of the couple ( ḟ , f ) versus x N within the Helmholtz ensemble.
A discussion concerning the duality and some conclusions on possible perspectives (Section 6) close the paper.

Configurational Partition Functions and Force-Extension Relations in Gibbs and Helmholtz Ensembles
The purpose of the present Section is to introduce a quite simple model which has the advantage to be analytically solvable for both the Gibbs isotensional ensemble and the Helmholtz isometric ensemble.The related mathematical analysis yields closed form expressions, which are beneficial to the thorough understanding of the physics of bistability (or, more generally, multistability) in complex systems, such as macromolecules of biological origin.
We consider a one-dimensional system composed of N elements with mechanical bistability, connected in series to compose a chain.Each element of the chain is therefore represented by a symmetric potential energy function U(x) showing two minima at x = ± (see Figure 1).As already described in the Introduction, to perform an analysis of the system reduced to essentials, we introduce a discrete variable y, which behaves as a spin, in place of considering the original bistable potential function represented in Figure 1 (blue dashed line).This spin variable pertains to the phase space of the system and, therefore, is a standard variable of the equilibrium statistical mechanics.The variable y assumes its values in the set S = {±1} and is used to identify the basin or well explored by the system.In conclusion, the original bistable energy function is substituted with the simpler mathematical expression The potential energy in Equation ( 1), by varying the value of the spin variable in S, generates the two parabolic wells represented in Figure 1 (red solid lines).While without an applied stretching the units are in each basin with the same probability (the average value of the end-to-end distance is zero), an applied stretching induces a preferential direction in the extension of the chain.This stretching can be applied by imposing a force f (positive or negative) or prescribing the position x N of the last element of the chain.Of course, in both cases, the first element is always tethered at the origin of the x-axis.These two possible mechanisms of stretching generate different stochastic mechanical behaviors of the system, which can be studied by calculating the corresponding configurational partition functions.

−`+0
x Bistable symmetric potential energy of a single domain (blue dashed line) and its approximation by means of four parabolic profiles (red solid lines).

The Gibbs Ensemble
In this case, we apply the force f to the last unit identified by its position x N .The total potential energy of the system under the Gibbs condition (isotensional ensemble) is therefore given by where f is the force applied to the last element, x = (x 1 , x 2 , ..., x N ) (continuous variables) and y = (y 1 , y 2 , ..., y N ) (discrete variables).For this system, we can define the configurational partition function Z G , as follows where the variable x is integrated whereas y is summed.We can now substitute Equation (2) in Equation (3).To evaluate the integral we apply the change of variables The change of variables within the multiple integral is implemented here by simply recalling that d x = Jd ξ.In this expression, the quantity J is the so-called Jacobian of the transformation defined as is the matrix of the first order partial derivatives ∂x i /∂ξ j .It can be easily proved that J = 1 for the proposed change of variables and, therefore, we finally get d ξ = d x, which strongly simplifies the calculation.Hence, we get where the integral I(y, f ) is defined as and it can be calculated in closed form by means of the well-known expression We eventually obtain the result Coming back to the configurational partition function, we have It is important to remark that within the Gibbs ensemble the elements of the chain do not interact and this point leads to a configurational partition function which is in the form of a power with exponent N.
The extension of the chain can be directly calculated through the expression x N = −∂U tot /∂ f and its average value is therefore x N = −∂U tot /∂ f .It can be simply evaluated by means of the configurational partition function, as We can also calculate the average value of the spin variable y = y i ∀i, which is independent of the element considered in the chain and is given by By combining Equation ( 9) with Equation ( 10), we immediately obtain This constitutive equation represents a spring-like behavior with an equilibrium length directly modulated by the average value of the spin variables.
An application of Equations ( 9) and ( 10) can be found in Figure 2. The force-extension curves have been plotted with dimensionless quantities and only one parameter defines the shape of the response, namely the elastic constant taken here into consideration through the dimensionless ratio k 2 k B T .It represents the ratio between the elastic (enthalpic) energy and the thermal energy.In these force-extension curves, we note a force plateau (for f = 0) corresponding to the synchronized switching (sometimes called cooperative) of the N units.This behavior is confirmed by the average spin variable (which is independent of k 2 k B T ), showing a transition from −1 to +1, at the same threshold force f = 0 as the previously mentioned plateau.This force plateau is the classical result of force-spectroscopy experiments conducted with soft devices [62][63][64][65][66][67].

The Helmholtz Ensemble
We can now introduce the second boundary condition corresponding to the Helmholtz ensemble.For imposing the isometric conditions, we consider the chain of bistable units with the two extremities tethered at the points x 0 = 0 and x N = x, respectively.The total potential energy of the system can be therefore written as where x N = x is the fixed extremity of the chain, x = (x 1 , x 2 , ..., x N−1 ) (continuous variables) and y = (y 1 , y 2 , ..., y N ) (discrete variables).In Equation ( 12) the potential energy U(x, y) of a single element is given in Equation ( 1).The configurational partition function of this system can be written as It it now important to remark that the isometric condition x N = x impedes the direct evaluation of the integral in Equation ( 13), which becomes considerably difficult.The solution to this problem can be found by drawing a comparison between Equations ( 3) and ( 13), eventually leading to the following important property: The two configurational partition functions Z G and Z H are related through a bilateral Laplace transform, as follows Moreover, if we let f = −iωk B T, we simply obtain which can be interpreted by affirming that the Fourier transform of Z H gives the analytical continuation of Z G on the imaginary axis of the complex plane.Exploiting this point, we can directly invert the Fourier transform, eventually obtaining Interestingly enough, we proved that the response of the system under the Helmholtz isometric ensemble can be analyzed through Equation (16), which considers as a starting point, the configurational partition function of the Gibbs isotensional ensemble.Anyway, from Equation (8), we have By using the Newton's Binomial Theorem we obtain from Equation ( 17) the partial result To go further, the integral in Equation ( 19) can be done with the help of the standard expression eventually obtaining It is interesting to observe that the isometric configurational partition function here obtained can not be stated in power form (with exponent N).This point suggests that under the Helmholtz ensemble there is an effective interaction among the elements of the chain.The origin of this interaction is not explicitly defined in the potential energy of the system (i.e., in the bistable character of the units), but comes from the specific boundary conditions characterizing the Helmholtz ensemble.Indeed, the isometric conditions fix the end-to-end distance by generating an effective interaction among the extensions of the units.Now, we can evaluate the average value of the overall force f = −k B T∂/∂x(log Z H ) applied to the system and the average value of the spin variables y = 1 N ∑ N i=1 y i describing the transitions, as follows and An example of application of Equations ( 22) and ( 23) can be found in Figure 3, where we show the force extension response and the average spin variable for the Helmholtz ensemble.As before, the force-extension curves have been plotted with dimensionless quantities and only one parameter defines the shape of the response, namely the elastic constant taken here into consideration through the dimensionless ratio k 2 k B T .We observe that the force extension curve is composed of a number of peaks corresponding to the non-synchronized (sequential) switching of the units.Sometimes, this behavior is called non-cooperative in order to underline the independent transitions of the units.This is confirmed by the step-wise curve representing the average spin variable versus the chain extension.Each step corresponds to the switching of a unit induced by the increasing extension of the chain.This behavior agrees with previous theoretical and experimental results obtained with hard devices [66][67][68][69][70][71][72][73].

Complete Probability Densities in the Gibbs and Helmholtz Ensembles
The results found in the previous Section concerning the Gibbs and Helmholtz partition functions and mechanical-configurational responses have been discussed for different systems in the scientific literature concerning the thermodynamics of bistability and the folding-unfolding processes.As we will show below, they represent the basis for investigating the behavior of these systems in more detail.In particular, we are interested not only in the average value of the fluctuating quantities, but also in the comprehensive statistical behavior described by the complete probability densities.The knowledge of these more refined quantities allows for the determination of expected values of higher order such as variances, covariances and so on, very important to fully characterize the statistical properties of these systems.We define here the probability density of the system in the whole phase space within both the Gibbs and the Helmholtz statistical ensembles.These results will be used in the following sections to find the probability density of the specific quantities characterizing the Gibbs and Helmholtz statistical ensembles.
Concerning the Gibbs ensemble, we can define the total energy of the system as where v i is the velocity of the i-th particle of the chain and v, x, y ∈ N while f ∈ .The complete probability density in the phase space is therefore given by the canonical distribution where the term 2πk B T m N has been added to normalize the kinetic part of the Boltzmann factor and the configurational partition function Z G ( f ) is given in Equation (8).Of course, we have that Similarly, for the Helmholtz ensemble we can define the total energy as where, as before, v i is the velocity of the i-th particle of the chain and v, x ∈ N−1 , y ∈ N while x N ∈ .In this case, the complete probability density in the phase space is given by the canonical distribution where the term 2πk B T m N−1 has been added to normalize the kinetic part of the Boltzmann factor and the configurational partition function Z H (x N ) is given in Equation (21).Of course, we have that The two probability densities here described will be used to obtain a full statistics representing the behavior of the two isotensional and isometric ensembles.

Probability Density of the Couple ( ẋN , x N ) versus f within the Gibbs Ensemble
Since the force f is imposed within the Gibbs ensemble, we can measure the extension of the chain which is a random variable that must be defined by its probability density in order to have a complete description of its behavior.Here, for the sake of completeness, we elaborate the probability density G ( ẋN , x N ; f ) for the couple ( ẋN , x N ), where we defined ẋN = v N .In this case, to obtain this probability density we have to sum or to integrate all the variables different from v N and x N in the complete density probability defined in Equation (25).It means that we can write Now, it is not difficult to recognize that the integral over the positions x 1 , ..., x N−1 immediately leads to the configurational partition function of the Helmholtz ensemble and the integral over the velocities v 1 , ..., v N−1 can be directly calculated with the classical Gaussian integral.Eventually, we obtain This is the most important result of this section and represents the probability density of the couple ( ẋN , x N ) for any value of the applied force f within the Gibbs ensemble.
We remark that this probability density can be factorized in two terms representing the density of ẋN and the density of x N .The first factor simply corresponds to the Maxwell distribution for the one-dimensional velocity On the other hand, it is interesting to observe that the second configurational term depends on the ratio between the two partition functions This configurational density is correctly normalized because of the Laplace integral relationship between Gibbs and Helmholtz partition functions, reported in Equation ( 14).The explicit form of G ( ẋN , x N ; f ) can be found by using the results given in Equations ( 8) and (21).The substitution yields the final expression By means of this expression, we can give another proof of the result giving the average value of x N .Indeed, we can write Now, the first integral is equal to 1 and the second one can be elaborated as follows By using again the Laplace integral relation between Gibbs and Helmholtz partition functions, reported in Equation ( 14), we easily get which is the well-known thermodynamic result.An example of application of the results obtained in the present Section is given in Figures 4-7.Since the kinetic component G ( ẋN ) is simply given by the Maxwell distribution, we focus our attention to the configurational part given by G (x N ; f ).Accordingly, in Figures 4 and 5, we show a three-dimensional and a two-dimensional representation of the Gibbs density as function of x N for a given applied force f .These results are represented for four different levels of thermal agitation in order to understand the effects of the disorder on the switching behavior between the states.The parameters used in this study are N = 5, = 1 (a.u.), k = 15 (a.u.) and k B T = 0.7, 1.4, 2.1, 2.8 (a.u.).It is interesting to observe that, in spite of the simple shape of the force-extension response characterized by a force plateau at f = 0 with a synchronized switching of the units, the probability density of the quantity x N is multimodal for the force range characterizing the transition region.Indeed, in order to obtain the probability density of x N for a given applied f we have to section the plots in Figures 4 and 5 with a plane parallel to the x N -axis and, at the same time, perpendicular to the f -axis.So doing, in the central transition region, we can observe the emergence of a series of peaks in the probability density confirming its multimodal character.This can be observed in Figure 6, where we plotted several curves G (x N ; f ) (see Equation ( 33)), for different values of the applied force f .We can observe the symmetric and multimodal profile of the probability density for f = 0 (at the center of the transition region) and the asymmetric and monomodal shape of the density for a large applied force (out of the transition region).We remark the multimodal character of the probability density of x N in spite of the simple force plateau observed in the force-extension response.To conclude this analysis, we underline that the knowledge of the full statistics for the system allows us to determine all possible expected values.As an example, we show in Figure 7 the behavior of the variance of the position in terms of the applied force f and the thermal energy k B T. We note that the variance is higher in the transition region, where the two states of each unit can coexist.Moreover, we observe a larger variance for higher temperatures, as expected.Finally, we also note that the multimodal character of the probability density is smeared out by the integration process applied to calculate the variances.This behavior will be shown to be dual with respect to the Helmholtz ensemble response, which is the subject of the next section.

Probability Density of the Couple ( ḟ , f ) versus x N within the Helmholtz Ensemble
The problem of finding the probability density for f and ḟ when x N is imposed is more complicated since, in this case, the variables f and ḟ do not belong to the phase space and, therefore, we can not integrate the superfluous variables in order to get the searched density.To cope with this problem, we first introduce the standard technique to deal with a function of random variable.We suppose to have two random variables x and y, linked by a function y = g(x).If F x (x) and f x (x) are distribution function and probability density of the random variable x, we search for the same quantities F y (y) and f y (y) for y = g(x).We use the symbol ξ for the elements of the probability space and we can write and Moreover, we can state that where 1(z) is the Heaviside step function.Therefore, we can obtain the probability density of y = g(x) by differentiation where we have introduced the Dirac delta function δ(z).This method based on the delta functions can be used to approach the problem of finding the Helmholtz probability density.To apply this technique, we need to write the variables f and ḟ in terms of the variables of the phase space of the system.Given the total potential energy Now, given the complete probability density ρ H ( v, x, y; x N ), we can obtain the density for the desired variables f and ḟ as follows This expression can be simplified delivering where we used the notation Z H (x N ) = Z N H (x N ) in order to specify that the Helmholtz partition function corresponds to a system with N units.Indeed, in the following calculations, we will also need the same function calculated for a system with N − 1 units.The elaboration of H ( ḟ , f ; x N ) can be continued as follows where we used the property δ(ax) = 1 |a| δ(x).We remember now that +∞ −∞ exp(−αx 2 ) = π α for α > 0, we perform the integrals of the delta functions over v N−1 and x N−1 , and we get We can now recall the explicit definition of Z N H (x N ) (see Equation ( 13)), and we also introduce the exact expression for So, in Equation ( 48), we can identify the partition function or, equivalently, This is the final result for the probability density within the Helmholtz ensemble.It is interesting to observe that it can be written in terms of the two partition functions Z N H (x N ) and Z N−1 H (x N−1 ), corresponding to systems of size N and N − 1, respectively.
We can split this probability density in two independent components describing separately ḟ and f , as follows and we can prove the normalization of the two results.For the first density function H ( ḟ ), the normalization directly comes from the classical integral +∞ −∞ exp(−αx 2 ) = π α for α > 0. For proving the normalization of the second density function H ( f ; x N ), we have to study the integral To do this, we observe that from Equations ( 49) and ( 50) we easily get the relation which can be also written as We can then make the changes of variables η a + = ξ and η b − = ξ, leading to Now, by letting x − ξ = f /k we eventually obtain that Finally, this result proves that H ( f ; x N ) is correctly normalized, being true that We also prove that we can re-obtain the well known expression for the average value of the force in the Helmholtz ensemble.To do this we consider the expression and we apply the change of variable x N − f /k = ξ leading to where we used Equation (59).As a conclusion, we proved that the classical thermodynamic relation for the average value of the force is consistent with our development.
A numerical application of the results proved in this section can be found in  Similarly to the Gibbs analysis, also in this case, we observe that the kinetic part of the probability density H ( ḟ ) is a simple Gaussian function and therefore we study in more detail the configurational density H ( f ; x N ).Coherently with this planning, in Figures 8 and 9, we show the three-dimensional and the two-dimensional representation of the Helmholtz density as function of f and for a prescribed extension x N .As before, the results have been obtained for four different temperatures to observe the effects of the thermal agitation on the transition processes.The parameters used in this study are the same already adopted for the Gibbs analysis, namely N = 5, = 1 (a.u.), k = 15 (a.u.) and k B T = 0.7, 1.4, 2.1, 2.8 (a.u.).We give here a description of the behavior of the system within the Helmholtz ensemble which is exactly dual with respect to the response of the Gibbs ensemble.Indeed, we observe that in spite of the saw-tooth shape of the force-extension response, the probability density of f is quite always monomodal.More precisely, it can be bimodal only with some sets of parameters and only for forces being in the transition region between two peaks of the force-extension curve.Anyway, we can affirm that this density is monomodal in the most cases of practical interest.To better explain this point, we observe that in order to obtain the probability density of f for a prescribed x N , we have to section the plots in Figures 8 and 9 with a plane parallel to the f -axis and, at the same time, perpendicular to the x N -axis.By performing this operation, in spite of the complex shape of H ( f ; x N ), we get monomodal functions (with the exceptions discussed above).This can be observed in Figure 10, where we plotted several curves H ( f ; x N ) for different values of the prescribed extension x N .As before, we remark that the knowledge of the full probability density for the Helmholtz case can be used to determine the expected values of higher order.As an example, in Figure 11 we plotted the variance of the force f , necessary to impose the extension x N .Interestingly enough, the variance is an increasing function of the temperature, as expected, and shows some peaks in correspondence to the switching of state of each unit.This is coherent with the general idea that the variance of the physical quantities is larger in proximity to a phase transition.Again, we underline the dual behavior of the Gibbs and Helmholtz ensembles.Indeed, while the variance for the Gibbs case is given by a single peak corresponding to the synchronized transition of the units, for the Helmholtz ensemble we have a peak for each transition, underlying the sequential behavior of this process.

Discussion and Conclusions
In this work we considered the comparison of Gibbs (isotensional) and Helmholtz (isometric) ensembles of the (equilibrium) statistical mechanics in the context of the stretching of chains of bistable units.The thermodynamics of the force-extension relations leads to different responses of the two ensembles for small systems, i.e., far from the thermodynamic limit.In particular, the Gibbs response is characterized by a force plateau corresponding to the synchronized transitions of the units, whereas the Helmholtz response can be viewed as a saw-tooth curve representing the sequential transitions of the units.We remark that, when the number of units approaches infinity, the two ensembles become equivalent from the statistical point of view and, therefore, the two Gibbs and Helmholtz force-extension responses become coincident.This general picture, well known in the context of the thermodynamics of small systems, has been widely confirmed experimentally by means of the force spectroscopy methodologies (see Introduction).
From the theoretical point of view, this scenario has been complemented here by introducing a method to elaborate the full statistics of these processes at thermodynamic equilibrium, i.e., by the calculation of the probability density of the fluctuating quantities and not only of their average values.The added information is useful to draw full comparisons with experiments and to extract more statistical features valuable from the theoretical point of view.As an example, the knowledge of the complete probability density can be used to evaluate expected values of higher order such as variances, covariances and so forth.Concerning the comparison with experiments, the devices today available to observe the mechanical response of macromolecules (force spectroscopy tweezers) are very refined and allow not only for the measurement of the average values of the main quantities but also to probe the distributions of the same quantities.This can be done by collecting the information of many trajectories of the system and to extract from those data the statistical picture of the system evolution.It means that it is important to update the theoretical devices in order to be able to calculate the probability density of the fluctuating quantities in terms of the deterministic applied ones.
Within the Gibbs ensemble, we apply a deterministic force and we measure a stochastic extension.So, we developed here a method to give the probability density of this extension and its rate with respect to the time.On the other hand, within the Helmholtz ensemble, we prescribe a deterministic extension and we measure a stochastic force.Therefore, we obtained in this work the probability density of the force and its time derivative.It is interesting to observe that in both cases these probability densities can be always written in terms of a combination of the two Gibbs and Helmholtz partition functions.This is a typical outcome in statistical mechanics, where all relevant quantities are typically written by means of the partition functions.We remark that, in the case of the number of units approaching infinity, we have the ensemble equivalence as above said.It means that the force-extension curves are the same for both ensembles but the probability densities are not the same because are simply defined on different variables.
The results obtained for the specific case of a chain of bistable elements show the emergence of an intriguing duality between the two ensembles.For the isotensional condition, the force-extension curve is monotone with a characteristic force plateau and the density (x N ; f ) is multimodal in the transition region (near x N = 0 and f = 0).Conversely, for the isometric condition, the force-extension curve is composed of a series of peaks while the density ( f ; x N ) is simply monomodal.This duality is also reflected in the behavior of the variances of these processes.In the Gibbs ensemble we obtained a monomodal variance σ 2 x with the symmetric peak at f = 0, whereas in the Helmholtz ensemble we obtained a multimodal variance σ 2 f with a peak for each transition value of x N .Of course, the peaks of variance must be explained through large fluctuations characterizing the switching of the units states (classically, it is typical for the phase transitions).
To go further with this analysis, in the next future we will take into consideration the case of bistable elements with two potential wells at different equilibrium length (as considered in this paper) and different equilibrium energy (here we supposed the same energetic level for the two basins).The introduction of the energy difference ∆E between the states will be useful to describe more realistic systems, such as protein domains and other macromolecules of biological origin.Another perspective concerns the consideration of the out-of-equilibrium thermodynamics useful to evaluate the dynamics (the time evolution) of the introduced probability densities, with application to the interpretation of force spectroscopy experiments.To do this, we plan to use our spin variables coupled with a first order dynamics governed by the Kramers rates, which depend on the energy barrier between the wells.
To conclude, it is important to underline that the thermodynamics of small systems and bistable-multistable systems is relevant not only for the studies concerning macromolecules and biophysics but only for several applications to nanoscience and nanotechnology, namely for the better understanding of plasticity, hysteretic behaviors and martensitic transformations in solids, micro-and nano-magnetism, ferromagnetic alloys, nano-indented substrates, bistable nanosystems for energy harvesting and transport phenomena in bistable nano-systems such as, e.g., tunnel effect transistors.

Figure 2 .
Figure 2. Average force-extension curves and average spin variables (plotted by means of dimensionless quantities) for the Gibbs ensemble with N = 5 and k 2 k B T = 10, 15, 30, 100.

Figure 3 .
Figure 3. Average force-extension curves and average spin variables (plotted by means of dimensionless quantities) for the Helmholtz ensemble with N = 5 and k 2 k B T = 10, 15, 30, 100.

Figure 6 .Figure 7 .
Figure 6.Examples of multimodal curves obtained through the Gibbs density G (x N ; f ) (see Equation (33)).On the left panel, the two-dimensional representation of the Gibbs density is shown with the cuts corresponding to the curves plotted on the right panel.We used N = 5, = 1 (a.u.), k = 15 (a.u.), k B T = 1 (a.u.) and different values of the applied force f , as indicated in the legend.

Figure 10 .
Figure 10.Examples of monomodal curves obtained through the Helmholtz density H ( f ; x N ) (see Equation (54)).On the left panel, the two-dimensional representation of the Helmholtz density is shown with the cuts corresponding to the curves plotted on the right panel.We used N = 5, = 1 (a.u.), k = 15 (a.u.), k B T = 1 (a.u.) and different values of the prescribed position x N , as indicated in the legend.