Ergodic Algorithmic Model (EAM), with Water as Implicit Solvent, in Chemical, Biochemical, and Biological Processes

For many years, we have devoted our research to the study of the thermodynamic properties of hydrophobic hydration processes in water, and we have proposed the Ergodic Algorithmic Model (EAM) for maintaining the thermodynamic properties of any hydrophobic hydration reaction at a constant pressure from the experimental determination of an equilibrium constant (or other potential functions) as a function of temperature. The model has been successfully validated by the statistical analysis of the information elements provided by the EAM model for about fifty compounds. The binding functions are convoluted functions, RlnKeq = {f (1/T)* g(T)} and RTlnKeq = {f (T)* g(lnT)}, where the primary linear functions f (1/T) and f (T) are modified and transformed into parabolic curves by the secondary functions g(T) and g(lnT), respectively. Convoluted functions are consistent with biphasic dual-structure partition function, {DS-PF} = {M-PF} · {T-PF} · {ζw}, composed by ({M-PF} (Density Entropy), {T-PF}) (Intensity Entropy), and {ζw} (implicit solvent). In the present paper, after recalling the essential aspects of the model, we outline the importance of considering the solvent as “implicit” in chemical and biochemical reactions. Moreover, we compare the information obtained by computer simulations using the models till now proposed with “explicit” solvent, showing the mess of information lost without considering the experimental approach of the EAM model.


Introduction
We accept in advance the suggestion of Lambert [1,2], who led a campaign against the theory which considers entropy as a parameter of disorder of the system. According to Lambert, entropy is a parameter of energy dispersion. It can be considered a parameter of disorder in as much as disorder is associated with energy. For instance, the sand of a beach is completely disordered unless at zero entropy, because there is no energy. By accepting the proposal of Lambert, we agree to defining free energy as a sum of two terms, Intensity Entropy (−∆H/T) and Density Entropy (∆S). The former function indicates the dispersion of energy produced by increasing temperature, T, which increases the velocity of the molecules, while the latter function indicates the dispersion of energy by increasing the dilution of the chemical species. We analysed the thermodynamic properties of hydrophobic hydration processes in water in a set of articles [3][4][5][6][7][8][9] and proposed the Ergodic Algorithmic Model (EAM) for maintaining the thermodynamic properties of any hydrophobic hydration reaction at a constant pressure.
If we want to obtain the thermodynamic functions of a chemical process, we plot the experimental data of an equilibrium constant RlnK dual , or other potential function RlnK eq , measured at different temperatures, as the function of (1/T), obtaining the convoluted parabolic binding function α) RlnK dual = {f (1/T)* g(T)}. If we plot the values of RTlnK dual as the function of T we obtain the convoluted parabolic binding function β) RTlnK dual = {f (T)* g(lnT)}. The binding functions are convoluted functions, where the primary linear functions f(1/T)and f(T)are modified and transformed into parabolic curves by the secondary functions g(T) and g(lnT), respectively. The constant factor R (gas constant) is introduced as a normalisation constant, to guarantee that we are referring to an Avogadro number of particles. Convoluted functions are consistent with biphasic dual-structure partition function, {DS-PF} = {M-PF} · {T-PF} · {ζ w }, composed by ({M-PF} (Density Entropy), {T-PF}) (Intensity Entropy), and {ζ w } (Implicit Solvent). In contrast, the functions obtained by simulation methods do not consider the free energy as a thermodynamic potential function (−∆G/RT) = (−∆H/RT) + (∆S/R) composed by two terms (Intensity Entropy (−∆H/RT) and Density Entropy (∆S/R)), specific for each hydrophobic hydration process. Simulations, referred to a conventional mono-phasic potential function or pseudo-free energy function, as proposed by other researchers, ignore the existence of the different structures of Intensity Entropy (−∆H/RT) and Density Entropy (∆S/R). The formal functions used in simulations are inconsistent with the dual biphasic structure of every experimental hydrophobic hydration system, with partition function composed by the product of a partition function {M-PF} · {T-PF} multiplied by the partition function of the implicit solvent ({ζ w }). The list of thermodynamic information elements provided by the Ergodic Algorithmic Model (EAM) and calculated by us for about fifty compounds was analysed by statistical analysis methodologies and successfully validated [8].

Thermal Equivalent Dilution (TED): Ergodicity
The Ergodic Algorithmic Model (EAM) is a set of mathematical equations whereby it is possible to extract from a series of experimental determinations of equilibrium constants, lnK, or other equivalent potential functions measured at different temperatures, the thermodynamic properties of any hydrophobic hydration reaction, by exploiting the Thermal Equivalent Dilution (TED) principle [3][4][5][6][7][8][9]. The TED principle can be represented by the equality n w C p,w dlnT = Rdlnd id,A where d id,A = 1/x A is ideal dilution; x A is the molecular fraction of species A; C p,w = 75.36 J K −1 mol −1 is the molar heat capacity of liquid water and R = 8.314 J·mol −1 K −1 is the gas constant. The expression in Equation (1) can be considered as the ergodicity parameter of a system. With reference to Figure 1a, representing a series of hyperbolas: The EAM model consists of a correct determination of the coordinates of a point A (x,y), with X = exp(∆S A /R) and y = exp(−∆H A /RT). By referring one equation to a given process, we choose one hyperbole of the family, corresponding to a point P on the auxiliary axis z of the diagram. Then, by choosing ( Figure 1b) the value of ordinate y (exp(−∆H A /RT)) and the value of abscissa x (exp(∆S A /R)) of point A (or of point B, or C), we choose a point on that hyperbola, thus reaching all the information elements available. Instead, this point A (or B or C) cannot be chosen by computer simulation. This is the essential information element that we lose by employing computer simulations: we cannot fix a point referring to the specific system on the chosen hyperbola. On the auxiliary coplanar axis z, we read the scale of values of point P, exp(−∆Γ/RT) = (exp(−∆G/RT) cos 45 • , with x P = y P . By computer simulation, we might perhaps choose the correct hyperbola P, but we cannot fix any point A for a specific reaction, losing all the essential information elements carried by the coordinates of A, −∆H A /RT, Intensity Entropy and ∆S A /R, Density Entropy, respectively.

Thermodynamic Free Energy
The probability space of Figure 1, representing the exponential functions, generates the thermodynamic space of Figure 2 respectively, where the exponentials are represented in an orthogonal axis system. In this system, the abscissa x A represents Density Entropy and the ordinate y  This function assumes implicitly that (−∆H (P) /T) = (∆S (P) ) on the diagonal, thus cancelling any difference, between the enthalpy parameter (or Intensity Entropy) and the entropy parameter (or Density Entropy). The computer-simulated function, not considering the difference between the two terms, Intensity Entropy and Density Entropy, cancels all the information elements contained there, and characterises each specific reaction, A, B, or C, at points A, B, or C, respectively. By applying the Ergodic Algorithmic Model (EAM) to a set of potential parameters (e.g., RlnK dual or RTnK dual ), referring to one specific hydrophobic hydration process, experimentally determined at different temperatures, one can calculate the binding functions α) RlnK dual = {f(1/T)·g(T)} (4) and which are "convoluted" functions. By developing the convoluted functions, one can calculate the primary motive functions, f(1/T) and f(T), respectively, of the convolution: The two terms of each motive function are specific for each compound and suited to recover the information elements of each hydrophobic hydration process. Other additional important information elements are produced by analysing the secondary function g(T) and g(lnT) of each convolution, such as ±ξ w , the pseudo-stoichiometric number of water molecules W I . The pseudo-stoichiometric number ±ξ w is calculated from the curvature of the binding function, that is a concave function, with ∆C p,hydr = ξ w C p,w >0 for Class A, or a convex function, with ∆C p,hydr = -ξ w C p,w < 0 for Class B. ξ w = |n w | is the number of water clusters W I involved in each process, and C p,w = 75.36 J·K −1 mol −1 is the molar heat capacity of liquid water.
The difference between the functions calculated by the so-called "simulation" and the convoluted binding functions obtained by applying the Ergodic Algorithmic Model (EAM) can be expressed by stating that the function calculated by simulation is lnK simul = f (−∆G/T), without specifying coordinates of any point A, B, or C, etc. The assignment of these equations means specifying, in the hyperbolic diagrams, the coordinates of each point, A, B, or C, each point referring to one specific different hydrophobic reaction, with specific properties and specific interrelations between properties for each family of compounds.
The complete set of equations constituting the Ergodic Algorithmic Model (EAM) are reported in Appendix A. The information elements from EAM are shown in Table 1.
with C p,w = 75.36 J·K −1 mol −1 (molar heat capacity for liquid water) ∆C p,hydr = +ξ w C p,w (Class A, convex binding function) The simulated calculated function RlnK simul is not suited to the study of these systems because it refers to a model of a monophasic system, not conforming to the biphasic composition [7] with implicit solvent of diluted aqueous systems and of biological solutions. The inadequacy of the monophasic model for simulations leads to the loss of essential information elements and to a pseudo-free energy potential function.

Pseudo-Free Energy Function
C. Chipot and A. Pohorille (Eds.) [10] published a book with title: "Free Energy Calculations. Theory and Application in Chemistry and Biology", where they treat the problem of the calculation of free energy functions lnK simul = f (−∆G/T). They optimize the function exp(−∆G/RT), but forget that thermodynamic free energy should be written and calculated as exp(−∆G/RT) = exp(−∆H/RT) · exp(∆S/R) (8) with a distinction between the Density Entropy term (exp(∆S/R) and Intensity Entropy term (exp(−∆H/RT)). Therefore, the name "Free Energy Calculations" should be changed into that of "Conventional monophasic" partition function or into that of "pseudo-Free Energy" function. If one is searching for a correct free energy procedure, one should set an algorithm to calculate . This means that we can determine, in advance, the final value of the two separate integrals by processing by EAM a set of potential functions, experimentally determined at three, four or more temperatures, and, in force of the Thermal Equivalent Dilution (TED) principle, behaving as dependent upon dilution. Another point ignored by Chipot and Pohorille [10,11] is the introduction of the constant partition function ζ w = 1, referring to water as an implicit solvent at constant thermodynamic potential µ s . Following the EAM procedure means cutting off all the treatments of hydrophobic processes produced according to the methods suggested by Chipot and Pohorille [10,11].
Though NEW shares common theoretical foundations with FEP, the latter can be assumed as a limiting case of the former [22].
Unfortunately, none of the eminent researchers [12][13][14]16,18,[20][21][22][23] active in simulation calculations have considered the thermodynamic free energy as composed by two distinct terms, Density Entropy and Intensity Entropy. On the other hand, Chipot and Pohorille and none of their followers have produced various sets of experimental data satisfying the statistical analysis and presenting mean unitary values with extremely low standard deviation. This important and decisive result was presented by us [8]. Free energy in chemical thermodynamics is set in probability space as and in thermodynamic space as: The partition function of Equation (10) has a typical hyperbolic structure, wherein exp(−∆H/RT) is representative of Intensity Entropy and exp(∆S/R) is representative of Density Entropy (Figure 1), as different distinct properties of a thermodynamic system. By choosing this free energy partition function, we accept that the thermodynamic system has a specific structure, composed of two entropic terms, i.e., Entropy Intensity and Entropy Density, respectively. These systems are constituted by a binary structure composed by a solvent at constant thermodynamic potential µ s (implicit solvent) and a solution with ergodic properties. These systems are represented in probability space by a dual-structure partition function.

Explicit Solvent and Implicit Solvent
Chipot and Pohorille [10] assumed the model of "explicit" solvent for water for the calculation of the pseudo-free energy. The "explicit" solvent [24] for water is incompatible with the dual-structure partition function {DS-PF} = {M-PF} · {T-PF} · {ζ w } of the aqueous systems as proposed by EAM model [9]. {ζ w } represents the partition function of implicit solvent, at constant thermodynamic potential µ s . These words explain the difference between pseudo-free energy calculations and the Ergodic Algorithmic Model (EAM). Pseudo-free energy calculations are applied to "explicit" solvent molecular systems, with homogenous composition whereas EAM requires an "implicit" solvent model for a dual system with biphasic composition. The dual system is composed by solvent in excess ("implicit" solvent) and diluted solute, as proved by statistical analysis of errors, extended to a large population of hydrophobic hydration processes [8]. Specifically, we assume [9] that every aqueous solution, and particularly biological solutions, has a biphasic composition, constituted by a solvent at constant potential (i.e., "implicit" solvent) and a diluted solute. None of the systems studied by perturbation theory [10,11] and by free energy perturbation (FEP) calculation [15] with determination of gradients and integration thereof, i.e., the partition function {ζ w }, give any contribution to the free energy of the whole system, conforming to the characteristics of the "implicit" solvent. The assumption of the "explicit" solvent for these systems would mean introducing a specific variable partition function for the explicit solvent. This hypothetical function should be suited to calculate a variable free energy contribution by the solvent, i.e., by the component at constant potential, whose contribution to free energy is constantly zero. In other words, we should introduce a partition function for the explicit solvent with the properties of implicit solvent, which would mean accepting the implicit solvent as the correct choice.
Notwithstanding these drawbacks, the computer simulation method called FEP (Free Energy Perturbation) [19] seems to have reached excellent results, (see Schrödinger Platform). Starting from some basic assumptions, several eminent researchers have developed a series of computer programs whereby an ingenious framework of mathematical relationships has been built up. We are dubious that these computer simulation methods are reliable.

Free Energy: Intensity Entropy and Density Entropy Components
We endeavour now to discover whether it is possible to create a mathematical framework, alternative to the pseudo-free energy function, with structures compatible with a structure of free energy (−∆G/RT), considered as the sum of one Intensity Entropy term (−∆H/RT) and one Density Entropy term (∆S/R), with the "implicit" solvent model. The integrals d(−∆H/RT) and d(∆S/R)} of Equation (9) and RTlnK mot = (−∆H 0 ) + (∆S mot )T (see Equations (6) and (7)). Pseudo-free energy calculations are applied to "explicit" solvent molecular systems with homogenous composition, whereas EAM requires an "implicit" solvent model for a dual system, composed by solvent in excess (implicit solvent) and diluted solute. The authors of simulations focus on the application of free energy perturbation methods utilizing MD for the calculation of protein-ligand binding affinities in structure-based drug discovery projects. The simulations consider the basic expression or Zwanzig equation, where T is the temperature and k B is Boltzmann's constant. The triangular brackets < and > denote an average over a simulation run for state A. In practice, one runs a normal simulation for state A, but each time a new configuration is accepted, the energy for state B is also computed. The difference between states A and B may be in the atom types involved, in which case the ∆F obtained is for "mutating" one molecule onto another, or it may be a difference of geometry, in which case one obtains an enthalpy map along one or more reaction coordinates. In practice, the averages in Equation (14) are obtained using the so-called λ-dynamics: E(λ) = (1−λ) · E A + λ · E B , where λ changes smoothly from 0 to 1 (i.e., from state A to state B) in order to improve the sampling of states [25].

Procedure According to Ergodic Algorithmic Model (EAM)
The dual system of every hydrophobic hydration process is represented in EAM by a product of partition functions {DS-PF} = {M-PF} · {T-PF th } · {ζ w }, whereby the Motive Partition Function {M-PF}, referring to Density Entropy, is multiplied by a thermal partition function {T-PF}, referring to Intensity Entropy, and by {ζ w }, at constant potential (implicit solvent). In mathematical format, we consider a partition function with ζ w = 1. According to EAM, the integration of differential functions cannot be applied to free energy calculation because the chemical reacting component, represented by {DS-PF}, is ruled by binomial distribution, with the partition function composed by a limited sum of finite arithmetic elements. EAM, however, considers an infinitesimal statistical distribution of molecules within sublevels h j,i of each macrolevel H i (see Figure 3). It is possible to apply the integration to the changing of distribution of molecules among sublevels, ruled by the statistical distribution of a non-reacting subsystem. In fact, the passage from one level to another implies a redistribution of molecules among the levels H i to H i+1 of the transformed molecule. This level-to-level passage can be a continuous process suitable for integration. The whole reaction also involves (Figure 4) a change in the multiplicity of states and this transformation of moles can be expressed by finite numbers, suitable for mathematical finite transformation corresponding to the coefficients of chemical equation. We can consider the contributions by integration terms as producing Intensity Entropy whereas the contributions by multiplicity mathematical terms as producing Density Entropy, respectively. The question now is: how can we exploit the computer programs on the market to calculate the Ergodic Algorithmic Model?
Instead of Zwanzig Equation (16), we should employ the equation In the former case, the solvent is considered at constant potential, i.e., as implicit solvent, whereas in the latter case the solvent is considered as explicit solvent.
In the calculation of the simulation function, we can consider the passage from one level H 0 to level H 1 with a gradual change in molecular structure from one to another shape, with a change in wave function. This passage consists of a change in Intensity Entropy: the whole mole passage of Intensity Entropy corresponds to the experimental determination of −∆H/T in the not Hoff function or similarly, in EAM to the slope of RlnK mot = −∆H/T + ∆S. The multiplicity within each macrolevel corresponds to Density Entropy change and is measured in experimental thermodynamic space by the ∆S component and in EAM by the slope of RTlnK mot = −∆H + T ∆S. There is correspondence between the two components of the equation −∆G/T = −∆H/T + ∆S and the terms calculated from the experimental data by applying RlnK mot in EAM.
The main achievements obtained by the Ergodic Algorithmic Model can be found in references [7][8][9].

Thermodynamic Functions and Information Elements
The development of the Ergodic Algorithmic Model (EAM) applied to a series of experimental data can show the relationships between thermodynamic functions and information elements.
The information necessary to construct the motive functions can be obtained from the definition of the primary functions f (1/T) of the binding function RlnK dual = {f (1/T)·g(T)} and by the primary function f (T) of the binding function RTln dual = {f (T)·g(lnT)}. The Ergodic Algorithmic Model (EAM) consists of the calculation and processing of the convoluted binding functions, to obtain the motive functions RlnK mot and RTlnK mot ( Table 2). The two motive functions contain the correct information elements concerning the leading hydration process of the reaction, together with the thermodynamic properties of the chemical process associated with hydration. In carboxylic acids, for example, the reaction associated with hydration is protonation of the base.   The correct procedure indicated by EAM is that followed by Talhout et al. when studying the binding to trypsin of Hexa-benzamidimium chloride. These authors began by determining by calorimetry the equilibrium constant of Hexa-benzamidimium chloride binding to trypsin [26] (Class B), at different temperatures for a series of hydrophobically modified benzamidinium chloride; then, they checked the results of simulations with the experimental findings.

Dual-Structure Partition Function. Molecular Structures
The conclusion of these comparisons between simulation procedures and the Ergodic Algorithmic Model (EAM) procedures is that computer simulations, although popular in the physical chemistry community, must be implemented by considering the diphasic structure of diluted systems. The reason for this criticism is that the mathematical functions employed by computer simulations and pseudo-free energy calculations are derived for a type of statistical ensemble not adequate for hydrophobic hydration processes. The statistical ensemble assumed by Chipot and Pohorille [10] to calculate RlnK simul = f (1/T) is monocentric and with a Boltzmann distribution with "explicit solvent", whereas EAM assumes that hydrophobic systems are represented by a dual-structure partition function We can obtain the information elements of Class A:
The same information elements are confirmed by the secondary function g(lnT) of the convoluted binding function RTln dual = {f (T) g(lnT)}.
In Class B, where we obtain the reaction of iceberg reduction (Figure 7), we obtain the information elements (a) (reaction: B {-ξ w W II (iceberg, solute)} → ξ w W I (solvent). (b) Concavity (∆C p,hydr = ξ w C p,w < 0). (c) Iceberg size and water stoichiometry: ξ w . From the coincidence of the value of entropy at zero iceberg ∆S 0 (ξ w = 0) = -86.4J K −1 mol −1 , (see reference [5]) with the Trouton constant ∆H eb /T eb = ∆S evap = +86.9 ± 1.4 J K −1 mol −1 , we can infer that the passage of gas molecules from the free state ( Figure 8) to the trapping into the solvent implies a change in entropy (as change in kinetic energy) that is exactly the opposite of that for the passage from liquid state to vapor state (i.e., configuration entropy change for the condensation of gas is equal to the opposite of evaporation, -∆S evap ).

Structures of Thermodynamic Functions
At the end of the analysis of the structure of the thermodynamic functions of hydrophobic hydration processes, we can conclude: (1) The simulated free energy functions RlnK simul , as calculated by any simulation functions, either FEP, or TI, MV, MC, etc., and applied to hydrophobic hydration processes are inadequate to represent the properties of biphasic systems with implicit solvent. Many important thermodynamic functions are ignored by these calculations and all the important information elements carried by these functions are completely lost. (2) The hydrophobic hydration systems have a biphasic composition, constituted by a "non-reacting" molecule ensemble, at constant potential formed by the solvent and by a "reacting" mole ensemble formed by diluted solutes. The distribution of molecules over the sublevels h i,j of each macrolevel H i follows a Boltzmann distribution, whereas the distribution of moles of reactants over macrolevels H i follows a binomial distribution. The Ergodic Algorithmic Model (EAM) is adequate to the biphasic composition of aqueous solutions and suited to evaluate correctly the many essential information elements.
The thermodynamic potential of each hydrophobic hydration process can be experimentally determined at different temperatures and processed by the Ergodic Algorithmic The statistical analysis of the unitary (for ξ = ±1) functions, which resulted to be constant in a statistically significant population of data, calculated by processing the curved binding functions, shows how the distribution of errors as a normal population validates the Ergodic Algorithmic Model (EAM) as the correct dual-structure model for the description of the thermodynamic properties of hydrophobic hydration processes. The statistical validation of the EAM model, discussed in ref. [8], is also reported in Appendix B. Another experimental method suited to the determinations of the information elements of a thermodynamic system in aqueous solution is calorimetry. Freire [27,28] showed how high affinity and selectivity, two essential properties of drugs, can be achieved by the optimization of enthalpic and entropic contributions to binding (the so-called "Thermodynamic Signatures"). If the calorimetric binding experiments are performed at a different temperature, the complete set of information elements can be obtained by applying the EAM model. This information is, therefore, essential not only from a speculative point of view, but also for application purposes, such as, for instance, lead optimization.
In any case, the statistical inference of Appendix B shows how the Ergodic Algorithmic Model (EAM) has a universal validity for every hydrophobic hydration process [8] and every hydrophobic hydration process must present the same behaviour. The functions required by the EAM must necessarily exist, for statistical inference. The first step is the search for an equilibrium constant or any appropriate potential function which is variable with temperature. The next step consists of a verification that the system presents the necessary ergodic properties. With these elements, one can determine the convoluted binding functions and then build up the EAM model.

Conclusions
The use of the Ergodic Algorithmic Model (EAM) with water as implicit solvent at constant thermodynamic potential µ s , to calculate the thermodynamic properties of hydrophobic hydration systems, is highly recommended, instead of computer simulations with explicit solvent, the evaluation of only binding affinity not being enough either from a speculative or practical point of view. Every hydrophobic hydration process necessarily has the properties of an ergodic model, statistically validated. For a correct thermodynamic analysis of any hydrophobic hydration process in chemical, biochemical, and biological systems, these properties must be searched for.