Dynamic Stability and Thermodynamic Characterization in an Enzymatic Reaction at the Single Molecule Level
DDynamic Stability and Thermodynamic Characterization in anEnzymatic Reaction at the Single Molecule Level
Mois´es Santill´an ∗ Centro de Investigaci´on y Estudios Avanzados del IPN, Unidad Monterrey,Parque de Investigaci´on e Innovaci´on Tecnol´ogica, 66600 Apodaca NL, M ´EXICO
Abstract
In this work we study, at the single molecular level, the thermodynamic and dynamic character-istics of an enzymatic reaction comprising a rate limiting step. We investigate how the stability ofthe enzyme-state stationary probability distribution, the reaction velocity, and its efficiency of en-ergy conversion depend on the system parameters. We employ in this study a recently introducedformalism for performing a multiscale thermodynamic analysis in continuous-time discrete-statestochastic systems.
Keywords:
Irreversible Thermodynamics; Non-equilibrium Steady State; Reaction Velocity; Ef-ficiency; Relaxation Rate
PACS numbers: 87.14.ej, 87.15.ad, 87.19.Pp ∗ Electronic address: [email protected] a r X i v : . [ q - b i o . S C ] M a y . INTRODUCTION The vast majority of processes taking place inside a cell consist of, or at least involveto some extent, chemical reactions. Moreover, most, if not all of the chemical reactionstaking place inside a cell are catalyzed by enzymes. Therefore, a profound understanding ofenzymes’ performance is necessary to better comprehend the processes of life.Homeostasis, understood as the coordinated physiological processes which maintain mostof the steady states in an organism, is regarded as a landmark concept in biology [1]. It can befound, to some extent, in all living beings, and allows them to perform in optimal conditionsdespite ever-changing surroundings and inputs. From a dynamical standpoint, homeostasisimplies the existence of a stable steady state [2–4]. Thus, the quality of homeostasis canbe measured by the volume of the steady-state basin of attraction in phase space and/orthe relaxation time with which the system returns to the steady state after a perturbation.Having a large basin of attraction is important because it allows the system to come back tothe steady state even in the face of large deviations. On the other hand, a rapid relaxationtime permits the system to quickly recover an optimal state after it is perturbed.Recent studies on finite-time thermodynamic engines and heat pumps have shown thattheir stability and their thermodynamic performance are often governed by the same param-eters. It has been observed that the system stability usually weakens as its thermodynamicproperties improve [5–17]. Consequently, these parameters need to be tuned to achieve anoptimal trade-off between favorable thermodynamic and dynamic properties. Similar studieson the stretch-reflex regulatory pathway and on a simple Brownian motor have confirmedthese findings [18, 19], in agreement with the notion that good design principles are usuallyshared by both artificial and biological systems [2]. If these results are of general applica-bility to a wide range of intracellular energy-converting processes, it would mean that themaintenance of the cell homeostatic state entails an expenditure of energy, which has tobe taken into consideration to understand how organisms adapt to a constantly changingenvironment.In the present work we study, at the single molecular level, the thermodynamic anddynamic characteristics of an enzymatic reaction. For the sake of simplicity we advocate toenzymatic reactions comprising a rate limiting step. We investigate how the stability of theenzyme-state stationary probability distribution, the reaction velocity, and its efficiency of2nergy conversion depend on the system parameters. We employ in this study a recentlyintroduced formalism for performing a multiscale thermodynamic analysis in continuous-time discrete-state stochastic systems.
II. MODELING ENZYMATIC REACTIONS
A simple but comprehensive enough model for an enzymatic reaction consists in picturingthe enzyme as undergoing the following series of chemical reactions [20]: first, a free enzyme E binds the substrate S ; then, the bound substrate is converted into the product P to formthe enzyme-product complex E P ; and finally, the product is released leaving the enzymefree to catalyze another reaction. This process can be summarized as follows using theconventional notation for chemical reactions: E + S (cid:10) E S (cid:10) E P (cid:10) E + P. Of all these reactions, the conversion of the substrate S into the product P ( E S (cid:10) E P )is in many cases the rate limiting process. Taking this into account and assuming thatthe substrate ([ S ]) and product ([ P ]) concentrations remain constant along the catalyticreaction, we can visualize a single enzyme as going through a series of transitions thatchange the enzyme state cyclically during the catalytic process; see Figure 1. There, theenzyme state is represented as ( i, j ), where index i denotes the i th cycle—that in which the i th product molecule is synthesized—while j = 1 , , E P , E + P , and E S . The assumption that [ S ] and [ P ] remain constant allows us to regard k + j and k − j ( j = 1 , ,
3) as pseudo first order reaction rates. Finally, the assumption that E S (cid:10) E P is the rate limiting process implies that vertical transitions in the scheme of Figure 1 aremuch faster processes than those involving changes in index i . III. PROBABILISTIC DESCRIPTION AND TIME-SCALE SEPARATION
Let us introduce a probabilistic description for an enzyme reaction at the single moleculelevel. And let us follow as well the approach in [21, 22] to simplify the model, takingadvantage of the assumed separation of time scales. Let P ( i, j ; t ) denote the probabilitythat the enzyme is in state ( i, j ) at time t . From the scheme in Figure 1, the chemical3 -1,2i-1,1i-1,3 i,2i,1i,3 i-1,2i-1,1i-1,3k +1 k -1 k +2 k -2 k +1 k -1 k +2 k -2 k +1 k -1 k +2 k -2 k +3 k -3 k +3 k -3 k +3 k -3 k +3 k -3 FIG. 1: Schematic representation of the various states available for an enzyme molecule and thetransitions rates among them, while the enzyme is catalyzing the synthesis of molecules P . Thestates are denoted as ( i, j ), where j = 1 , , E P , E + P , and E S , while the index i denotes the i th enzyme cycle: that in which the i th product molecule issynthesized and released. master equation for P ( i, j ; t ) consists of the following set of coupled differential equations: dP ( i, t ) dt = k +3 P ( i − , t ) + k − P ( i, t ) − ( k − + k +1 ) P ( i, t ) , (1) dP ( i, t ) dt = k +1 P ( i, t ) + k − P ( i, t ) − ( k − + k +2 ) P ( i, t ) , (2) dP ( i, t ) dt = k +2 P ( i, t ) + k − P ( i + 1 , t ) − ( k − + k +3 ) P ( i, t ) . (3)The probability that the enzyme is in state ( i, · ) at time t is given by P ( i ; t ) = (cid:88) j =1 P ( i, j ; t ) . (4)On the other hand, it follows from the definition of conditional probability that P ( i, j ; t ) = P ( j | i ; t ) P ( i ; t ) . (5)Add Equations (1)-(3) and use (4) and (5) to obtain dP ( i ; t ) dt = (cid:0) k +3 P (3 | i − t ) (cid:1) P ( i − t ) + (cid:0) k − P (1 | i + 1; t ) (cid:1) P ( i + 1; t ) − (cid:0) k +3 P (3 | i ; t ) + k − P (1 | i ; t ) (cid:1) P ( i ; t ) . (6)Differentiate (5) and assume a time-scale separation so that the transitions between states( i − ,
1) and ( i,
3) are much slower than all the other (that is, they are the rate limiting4teps along the reaction chain). Then, dP (1 | i ; t ) dt = k − P (2 | i ; t ) + k +1 P (1 | i ; t ) , (7) dP (2 | i ; t ) dt = k +1 P (1 | i ; t ) − k − P (2 | i ; t ) + k − P (3 | i ; t ) − k +2 P (2 | i ; t ) , (8) dP (3 | i ; t ) dt = k +2 P (2 | i ; t ) − k − P (3 | i ; t ) . (9)If we invoke once more the separation of time scales to assume that the fast dynamicsrapidly reach an equilibrium distribution while slow dynamics have not changed noticeably: dP ( j | i ; t ) /dt ≈ j = 1 , ,
3, we have that P (1 | i ; t ) (cid:39) P ∗ (1 | i ) ,P (2 | i ; t ) (cid:39) P ∗ (2 | i ) = K P ∗ (1 | i ) ,P (2 | i ; t ) (cid:39) P ∗ (3 | i ) = K P ∗ (2 | i ) = K K P ∗ (1 | i ) , where K j = k + j /k − j , j = 1 , ,
3. Finally, the normalization condition ( (cid:80) j =1 P ∗ ( j | i ) = 1)implies that P ∗ (1 | i ) = 11 + K + K K , (10) P ∗ (2 | i ) = K K + K K , (11) P ∗ (3 | i ) = K K K + K K . (12)Notice that P ∗ ( j | i ) ( j = 1 , ,
3) are all independent of i . Furthermore, it is straightforwardto prove that this stationary distribution satisfies the following relations: k +1 P ∗ (1 | i ) = k − P ∗ (2 | i ) , and k +2 P ∗ (2 | i ) = k − P ∗ (3 | i ) . This last result implies that the fast dynamics are in equilibrium only if each one of theunderlying chemical reactions is in equilibrium itself. Finally, substitution of (10)-(12) into(6) allows us to conclude that, when time-scale separation is possible and the enzyme statesare grouped as sketched in Figure 1, the system dynamics is that of a biased one-dimensionalrandom walk: dP ( x ; t ) dt = k + P ( i − t ) + k − P ( i + 1; t ) − (cid:0) k + + k − (cid:1) P ( i ; t ) , (13)where k + = γ K K K K + K K , and k − = γ
11 + K + K K , (14)5hile γ = k − and K = k +3 /k − .From (13), the slow-dynamics stationary distribution obeys P ∗ ( i ) = P ∗ ≡ constant . Observe that this distribution does not fulfill detailed balance since k + P ∗ ( i − − k − P ∗ ( i ) (cid:54) =0, unless k ∗ = k − . IV. THERMODYNAMIC STATE VARIABLES AND RELAXATION TO THESTATIONARY STATE
Following [22–24], the enzyme internal energy, entropy, and Helmholtz free energy canbe respectively defined as follows: U = − k B T (cid:88) i,j P ( i, j ; t ) log P ∗ ( i, j ) ,S = − k B (cid:88) i,j P ( i, j ; t ) log P ( i, j ; t ) ,F = U − T S = k B T (cid:88) i,j P ( i, j ; t ) log P ( i, j ; t ) P ∗ ( i, j ) . By substituting (5) into the above equations they can be rewritten as U = − k B T (cid:88) i P ( i ; t ) log P ∗ ( i ) − k B T (cid:88) i P ( i ; t ) (cid:88) j P ( j | i ; t ) log P ∗ ( j | i ) , (15) S = − k B (cid:88) i P ( i ; t ) log P ( i ; t ) − k B (cid:88) i P ( i ; t ) (cid:88) j P ( j | i ; t ) log P ( j | i ; t ) , (16) F = k B T (cid:88) i P ( i ; t ) log P ( i ; t ) P ∗ ( i ) + k B T (cid:88) i P ( i ; t ) (cid:88) j P ( j | i ; t ) log P ( j | i ; t ) P ∗ ( j | i ) . (17)Observe that this way of writing the thermodynamic state variables renders a natural sep-aration of contributions from the fast and slow dynamics.The first term in the right hand side of Equation (17) is nothing else but the Kullback-Leibler divergence between distributions P ( i ; t ) and P ∗ ( i ) and so it is positive defined andonly equals zero when the two distributions are identical. Similarly, the sum over j in thesecond term is the Kullback-Leibler divergence between P ( j | i ; t ) and P ∗ ( j | i ), it is positivedefined, and only equals zero when P ( j | i ; t ) = P ∗ ( j | i ). From these considerations, the valueof F can be used as an indicator of how far the system is from the stationary distribution.6fter differentiating Equation (17) and imposing the separation of time scales we obtain dFdt = Q kh − T σ, (18)where Q hk is known as the housekeeping heat and is given by Q hk = k B T (cid:88) i (cid:0) P ( i ; t ) k + − P ( i + 1; t ) k − (cid:1) log P ∗ ( i ) k + P ∗ ( i + 1) k − + k B T (cid:88) i P ( i ; t ) (cid:88) j =1 (cid:0) P ( j | i ; t ) k + j − P ( j + 1 | i ; t ) k − j (cid:1) log P ∗ ( j | i ) k + j P ∗ ( j + 1 | i ) k − j , (19)while σ is the entropy production rate: σ = k B (cid:88) i (cid:0) P ( i ; t ) k + − P ( i + 1; t ) k − (cid:1) log P ( i ; t ) k + P ( i + 1; t ) k − + k B (cid:88) i P ( i ; t ) (cid:88) j =1 (cid:0) P ( j | i ; t ) k + j − P ( j + 1 | i ; t ) k − j (cid:1) log P ( j | i ; t ) k + j P ( j + 1 | i ; t ) k − j . (20)Observe that both Q hk and σ have contributions from the slow (first term on the right handside) and fast dynamics (second term). However, the fast dynamics contribution to Q hk vanishes because P ∗ ( j | i ) complies with detailed balance, and so Q hk = k B T (cid:88) i (cid:0) P ( i ; t ) k + − P ( i + 1; t ) k − (cid:1) log P ∗ ( i ) k + P ∗ ( i + 1) k − . (21)This result is in complete agreement with the interpretation of Q hk as the energy that hasto be pumped into the system to drive the stationary state out from equilibrium (detailedbalance).In the present formalism, the enzyme molecule is implicitly assumed to be in equilibriumwith a thermal bath. Concomitantly, the proper thermodynamic description is that of theHelmholtz free energy. As seen in Equation (17), F ≥ F can be used as a measure of how distant thesystem state is from the stationary one, as we have previously asserted. Furthermore, it isnot hard to prove from (18), (20), and (21) that dF/dt ≤
0, and that dF/dt = 0 only when P ( i ; t ) = P ∗ ( i ) and P ( j | i ; t ) = P ∗ ( j | i ). Therefore, dF/dt can be interpreted as the rate ofrelaxation to the stationary distribution.We can see from (18), (20), and (21) that dF/dt can be decomposed into contributionsfrom the slow and fast dymamics: dFdt = ˙ F slow + ˙ F fast , (22)7here ˙ F slow = k B (cid:88) i (cid:0) P ( i ; t ) k + − P ( i + 1; t ) k − (cid:1) log P ∗ ( i ) P ( j + 1 | i ; t ) P ∗ ( i + 1) P ( j | i ; t ) , (23)and ˙ F fast = k B (cid:88) i P ( i ; t ) (cid:88) j =1 (cid:0) P ( j | i ; t ) k + j − P ( j + 1 | i ; t ) k − j (cid:1) log P ( j | i ; t ) k + j P ( j + 1 | i ; t ) k − j . (24)The relaxation of the fast dynamics subspace to the corresponding quasi-stationarystate—given by the slow-dynamics probability distribution—is determined by ˙ F fast . Accord-ing to the assumed separation of time scales, fast dynamics relaxation takes place withoutany noticeable change in the slow dynamics probability distribution ( P ( i ; t )) and, in thelong run, the slow dynamics are the ones that govern the system relaxation to the steadystate. That is, the system relaxation rate is well approximated by ˙ F slow .Under the assumption that the probability distribution P ( i ; t ) is slightly different fromthe stationary distribution P ∗ : P ( i ; t ) = P ∗ + (cid:15) ( i ; t ), with (cid:15) ( i ; t ) (cid:28)
1. Then, it is straightfor-ward to see from (23) that, in the neighborhood of the stationary distribution, the systemrelaxation rate (here denoted by ξ ) is proportional to—see equation (14): ξ ∝ k + − k − = γ K K K −
11 + K + K K . (25)To better understand the result in Equation (25) let us analyze the significance of pa-rameters K i ( i = 1 , , E P (cid:10) E + P ( K ), E + S (cid:10) E S ( K ),and E S (cid:10) E P ( K ), respectively. On the other hand, a reaction’s association rate ( K A )is related to its free energy change (∆ G ) by K A = exp( − ∆ G/RT ), with R the ideal gasconstant and T the absolute temperature. These considerations allow us to visualize anenzymatic reaction as a process occurring along an energy profile like the one pictured inFigure 2. In such scheme, the enzyme states correspond to the local minima of the energyprofile, and the transition probabilities are determined by the height of the energy barriers.The global free energy change (∆ G T ) of an enzymatic reaction is given by ∆ G T = ∆ G +∆ G + ∆ G = − RT log( K K K ). Therefore, since the presence of an enzyme does notchange ∆ G T : K K K = exp (cid:18) − ∆ G T RT (cid:19) ≡ constant . This restriction further implies that only two of the three K i constants are independent.Without loss of generality we shall consider that K and K are determined by the nature of8 nzyme States E ne r g y P r o fi l e (i,1)E P (i,2)E+P (i,3)E S (i,4)E P G G G G Δ G > 0 Δ G + Δ G + Δ G < Δ G > 0 Δ G < 0 FIG. 2: A cartoon representation of the Gibbs free energy profile for an enzymatic reaction. Theminima correspond to the states the enzyme goes through, while the transition probabilities aredetermined by the height of the energy barriers. The individual free energy changes betweenadjacent states can be either positive or negative, but the global free energy change has to benegative in order for the reaction to proceed forward. the enzyme catalyzing the reaction, while K = exp( − ∆ G/RT ) /K K . With this, Equation(25) can be rewritten as ξ ∝ γ exp( − ∆ G T /RT ) −
11 + K + K K . (26)Note that the relaxation rate is a monotonic decreasing function of both K and K . If wefurther take into consideration that K , K >
0, it follows that the relaxation rate can beincreased by making K and K as close to zero as possible, and therefore by making ∆ G and ∆ G as positively large as possible. In particular, the maximum value of ξ is attainedwhen K = 0, regardless the value of K . Thus, in order to increase the value of ξ it is moreimportant to increase the value of ∆ G than that of ∆ G .On the other hand, we can see from Equation (26) that the relaxation rate is a monotonicdecreasing function of ∆ G T . That is, the more energetically favorable the global functionis, the more stable the stationary distribution becomes. Similarly, ξ is proportional to γ = − k . Therefore, the system stationary distribution is more strongly stable when thesubstrate-to-product conversion process is more rapid.9 . REACTION VELOCITY AND EFFICIENCY OF ENERGY CONVERSION As we have seen, by exploiting the time-scale separation to simplify the reaction scheme,an enzyme can be modeled as a biased one-dimensional random walk with forward andbackward transition probabilities k + and k − , respectively—see Equations (13) and (14).Therefore, by taking into consideration that the stationary probability distribution P ∗ isconstant, the reaction velocity can be defined as the forward minus the backward fluxes: ν = P ∗ ( k + − k − ) = γP ∗ exp( − ∆ G T /RT ) −
11 + K + K K . (27)On the other hand, if we consider that a certain amount of energy is consumed during eachforward step, and that this energy is waisted during backward steps, the system efficiencyin the stationary state can be define as η = 1 − k − k + = 1 − − ∆ G T /RT ) . (28)Observe that ν > G T , and that ν = 0 when ∆ G T = 0. Inother words ∆ G T < η is also adecreasing function of ∆ G T , and η = 0 when ∆ G T = 0. That is, both the reaction velocityand its efficiency can be increased by making ∆ G T more negative.Note from equations (26) and (28) that η ∝ ξ . Therefore, the discussion regarding thedependence of ξ on K and K applies as well to η . In particular, we want to emphasize thatthe reaction velocity can be increased by making K and K as close to zero as possible (andthus by making ∆ G and ∆ G as large as possible). However, varying K is more importantsince the maximum velocity can be achieved by setting K = 0, regardless the value of K .Interestingly, the reaction efficiency is independent of K and K . Thus, given that anenzyme does not alter the global free energy change of the reaction it catalyzes, this resultimplies that a reaction efficiency is the same regardless whether it is catalyzed or not. VI. CONCLUDING REMARKS
In this work we have investigated the dynamic stability, as well as velocity and efficiency,of an enzymatic reaction with a rate limiting step, at the single molecule level. For this,we followed the ideas in [23, 24], and used a recently developed formalism for performing10ultiscale thermodynamic analysis on discrete-state, continuos-time, Markovian stochasticprocesses [22]. Our results can be summarized as follows:1. The dynamic and thermodynamic characteristics associated to the stationary proba-bility distribution are completely determined by the the Gibbs free energy changes ofthe enzymatic reaction steps: ∆ G ( E P (cid:10) E + P ), ∆ G ( E + S (cid:10) E S ), and ∆ G ( E S (cid:10) E P ). Observe that the energies ∆ G i ( i = 1 , ,
3) are not all independent becase (cid:80) i =1 ∆ G i = ∆ G T , and ∆ G T (the global free energy change) is not modified by theenzyme.2. The stationary probability distribution is stable and the corresponding relaxation rate( ξ ) is directly proportional to the global reaction velocity ( ν ).3. Both ξ and ν are decreasing functions of ∆ G T , and ξ, ν = 0 when ∆ G T = 0. Thus,the global reaction accelerates and the stationary probability distribution turns morestrongly stable as ∆ G T is more negative.4. Both ξ and ν are increasing functions of ∆ G and ∆ G . The relaxation rate and thereaction velocity achieve their maximum value in the limit ∆ G → ∞ , regardless thevalue of ∆ G . Contrarily, ξ and ν increase as ∆ G increases and converge to a valuethat depends on ∆ G as ∆ G → ∞ .5. The efficiency ( η ) is a function of ∆ G T , but it is independent of ∆ G and ∆ G . Inparticular, η is a decreasing function of ∆ G T , and η = 0 when ∆ G T = 0. That is, thereaction efficiency increases as ∆ G T becomes more negative.Regarding the novelty of the above results, it is worth pointing out that, while someof them (like the one stating that ν decreases as ∆ G T increases) are well known, there areothers (like those regarding the relaxation rate ξ ) which are new to the best of our knowledge.On the other hand, when the separation of time scales is not possible, the system dynamicand thermodynamic characteristics will be determined by the rate constants k + i and k − i ( i = 1 , , ξ , ν , and η , shall not necessarily be given in terms of theequilibrium dissociation constants K i = k − i /k + i , as in Equations (26)-(28). The reason beingthat the chemical reactions taking place between the synthesis of one product molecule andthe next won’t necessarily be close to chemical equilibrium. Some of the formerly enumerated11esults may still be qualitatively true for enzymatic reactions without a rate limiting step,though. However, the study of such reactions is beyond the scope of the present work. Inthe forthcoming paragraphs we analyze the implications of the results listed above; in theunderstanding that, unless otherwise stated, these implications may not be exclusive forenzymatic reactions with a rate limiting step.Although in an enzymatic reaction the global free energy change is predetermined by thesubstrate, the products, and the nature of their reaction, the actual shape of the energyprofile along the reaction coordinate depends on the enzyme structure. Our results confirmthat although the efficiency is not affected by the shape of the energy profile, the reactionvelocity and the strength of the stationary distribution stability can be highly improved byproperly shaping this profile. This behavior is contrary to that observed in other systems(like thermal engines) in which variation in some parameters makes the velocity and theefficiency change in opposite directions. In other words, in this case no trade-off betweenefficiency and reaction velocity (and stability strength) need to be looked for, regardingthe energy profile. We believe this may by one of the reasons why it has been possiblefor evolution to drive the structure of enzymes so the corresponding reaction velocity isincreased by several orders of magnitude.Another feature worth noticing is the fact that the more energetically unfavorable reac-tions E P (cid:10) E + P and E + S (cid:10) E S are, the faster the global enzymatic reaction is. Thisbehavior can be understood by looking at Figure (2). We can see there that large, positive∆ G and ∆ G values imply a very negative ∆ G . This further means that the backwardreaction E S (cid:41) E P is much less probable than the forward reaction E S (cid:42) E P . That is,the strategy to accelerate the global enzymatic reaction seems to be making the rate lim-iting step almost unidirectional, even though this implies that, in the rapid processes, thebackward reactions are more probable than the corresponding forward reactions. Given theprominent role the rate limiting step takes in the above discussion, we believe it does notapply to other types of enzymatic reactions. Acknowledgments
This work was partially supported by Conacyt, M´exico, Grant: 55228. The authoris thankful to the International Centre for Theoretical Physics for its support during the12ealization of this work. [1] F. L. Holmes, Hist. Philos. Life Sci. , 3 (1986).[2] H. Kitano, Nat. Rev. Genet. , 826 (2004).[3] H. Kitano, Mol. Sys. Biol. , 137 (2007).[4] A. Lesne, Biol. Rev. Camb. Philos. Soc. , 509 (2008).[5] M. Santill´an, G. Maya, and F. Angulo-Brown, J. Phys. D: Appl. Phys. , 2068 (2001).[6] L. Guzman-Vargas, I. Reyes-Ramirez, and N. Sanchez, J. Phys. D: Appl. Phys. , 1282(2005).[7] J. Chimal-Eguia, M. Barranco-Jimenez, and F. Angulo-Brown, Open Syst. Inf. Dyn. , 43(2006).[8] R. P´aez-Hern´andez, F. Angulo-Brown, and M. Santill´an, J. Non-equilib. Thermodyn. , 173(2006).[9] Y. Huang, D. Sun, and Y. Kang, J. Appl. Phys. (2007).[10] J. C. Chimal-Eguia, I. Reyes-Ramirez, and L. Guzman-Vargas, Open Syst. Inf. Dyn. , 411(2007).[11] W. Nie, J. He, B. Yang, and X. Qian, Appl. Therm. Eng. , 699 (2008).[12] W. Nie, J. He, and X. Deng, Int. J. Therm. Sci. , 633 (2008).[13] Y. Huang and D. Sun, Int. J. Refrig. , 483 (2008).[14] Y. Huang and D. Sun, J. Non-equilib. Thermodyn. , 61 (2008).[15] Y. Huang and D. Sun, Appl. Therm. Eng. , 1443 (2008).[16] Y. Huang, D. Sun, and Y. Kang, Appl. Therm. Eng. , 358 (2009).[17] Y. Huang, Energy Convers. Manage. , 1444 (2009).[18] R. P´aez-Hern´andez and M. Santill´an, Physica A , 3574 (2008).[19] M. Santill´an and M. C. Mackey, Phys. Rev. E (2008).[20] H. Qian and E. L. Elson, Biophys Chem , 565 (2002).[21] E. S. Zeron and M. Santill´an, J Theor Biol , 377 (2010).[22] M. Santill´an and H. Qian, Phys. Rev. E In Press (2011), arXiv:1003.3513v6.[23] M. Esposito, U. Harbola, and S. Mukamel, Phys. Rev. E , 031132 (2007).[24] H. Ge and H. Qian, Phys. Rev. E , 051133 (2010)., 051133 (2010).