Recent Developments on the Mechanism and Kinetics of Esterification Reaction Promoted by Various Catalysts

Esters have played a significant role in daily living and chemical industry, such as plasticizers, fragrance, adhesive and lubricants (Joseph et al., 2005; Mbaraka & Shanks, 2006; Krause et al., 2009; Martínez et al., 2011). The vast majority of esters can be prepared using esterification reaction in the chemical engineering industry. Esterification has acquired further improvement from the engineering side; this mainly depends on the research of esterification kinetics. On the other hand, the need to control chemical reactions at the molecular level, which depends critically on the catalytic mechanism, is rapidly increasing (Salciccioli et al., 2011).

The chemical reaction can be expressed as Figure 1.   (Streitwieser et al., 1985;Fei & Zhao, 2009 In addition, Streitweiser considered that step (2) is the rate-determining stage. The total reaction is: The result of the reaction is that the rupture of carboxylic acid acyl-oxygen bond occurs, and the hydroxy of acid is replaced by alkoxy, which is a nucleophilic substitution process of carboxylic acids.
Esterifications of primary alcohols, secondary alcohols with carboxylic acids comply with this mechanism usually. According to this mechanism, the intermediate with a tetrahedral structure is more crowd than reactants, which makes the structure of carboxylic acid and alcohol have a significant effect on the easy of esterification.
The activity of different structures of the alcohol and carboxylic acid in esterification obeys the following order: Alcohol: CH 3 OH > RCH 2 OH > R 2 CHOH; Acid: CH 3 COOH > RCH 2 COOH > R 2 CHCOOH > R 3 CCOOH (Solomons, 1986;March, 1992;TianQuan 1992;Fei & Zhao, 2009).It shows that straight-chain structure is easier than branched-chain structure to esterification, and the more the branched-chain is, the lower the rate will be.
Based on the above mechanism proposed by Streitweiser, Rönnback, et al. developed a esterification kinetic model of carboxylic acid with methanol in the presence of hydrogen iodide though isothermal batch experiments at 30-60℃. The catalyst concentration varied from 0.05 to 10.0 wt%. Because the proton-donation step (1) as well as the subsequent steps (3)-(5) is assumed to be rapid, the simplified mechanism shown in Figure 3: where r is the total rate, r 2 is the rate of step(2), A is CH 3 C(OH) 2 + and E is CH 3 COOCH 3 .
The kinetic and equilibrium parameters included in Equation (5) are estimated from experimental data with regression analysis. Simulation of the model with the estimated parameters revealed that Equation (5) can predict the experimental trends in the acidcatalyzed esterification correctly.

Single-molecule reaction
The tertiary alcohol is prone to generate carbocation under acidic condition. Therefore, the mechanism of tert-esterification is different from that of the primary and secondary alcohol, it follows the single-molecule reaction process(shown in Figure 4):The tertiary alcohol combines with protons to generate protonated alcohol(1); The protonated alcohol is got rid of a molecule water to produce tert-carbocation(2), which is very stable; Then the electrophilic attacking is taken place between the tert-carbocation and oxygen atom of carboxylate, the protonated ester(3) is yielded; The loss of a proton from the protonated ester gives the product(4) (Bart et al., 1994).

Large-steric reaction
Aromatic acid has serious steric hindrance. If there are both ortho-methyls, just as 2.4.6trimethyl benzoic acid, the alcohol molecule is so difficult to access to carboxyl that the esterification cannot occur. However, if 2.4.6-trimethyl benzoic acid is dissolved in the 100% sulfuric-acid solution, acylilum ion will be formed as shown in Figure 5. Then the added alcohol with the acylilum ion produced ester. The reaction will conduct smoothly (Streitwieser et al., 1985).
This kind of esterification could occur for the reason that carbon atom of the acylilum ion is sp 2 hybridized and coplanar with the benzene ring, Then alcohol molecule can be virtually unhindered to attack the acylilum ion from above or below of the molecular plane. Esterifications conducts with this mechanism are merely few.
The alkylation reaction of sulfuric acid with alcohol is irreversible (Aranda et al., 2008).Therefore the reaction rate can be expressed as: where k is rate constant and K is equilibrium constant; C A , C B , C W , C E is the concentration of carboxylic acid, alcohol, water and ester, respectively. In 2009, Fei, et al. studied on the esterification behaviors of neo-polyhydric alcohols with fatty acids catalyzed by sulfuric acid. It was found that the esterification follows the above mechanism.

Lewis acids
Lewis acid is a compound or ionic species that can accept an electron pair from a donor compound to form a Lewis adduct, such as ZnCl 2 , Mg (ClO 4 ) 2 , Sn[N(SO 2 -n-C 8 F 17 ) 2 ] 4 . It is potentially useful on the production of esters as a catalyst, the reason is that it can be easily separated from the reaction media (Sanchez et al., 1992;Cardoso et al., 2009); moreover, comparing with Bronsted acid, template effects are to be expected as Lewis acid sterically bulkier than a proton. In 1964, Anantakrishnan, et al. investigated the esterification reaction of Ac 2 O with 50% MeOH, EtOH, PrOH, and iso-PrOH, in 50% Me 2 CO or Dioxane as solvent, using Lewis acid (ZnCl 2 ) as catalyst.
Along with the scientific and technological progress and social development, more and more researchers have concerned about this issue, including the kinds of catalysts, the mechanism and kinetics of reaction. Thus, the contents containing the reaction mechanism and esterification kinetics have been significantly reviewed according to the type of catalysts as follows.

Metal chloride
Ethyl oleate was synthesized by the esterification of and ethanol catalyzed by SnCl 2 •2H 2 O . Under the circumstance of excess ethanol, the effects of the concentration of the catalyst and oleic acid, and temperature on the reaction rate were investigated. A related esterification mechanism was presented and described as follows: in presence of Sn 2+ ( SnCl 2 •2H 2 O) catalyst, the carbonyl of the fatty acid is polarized to activate of substrate, which makes the nucleophilic attack to the molecules by ethanol become more favorable. Cardoso et al. investigated the effect of different carbonic chain of alcohol (methyl alcohol, ethyl alcohol, n-propyl alcohol, n-butyl alcohol) on the conversion of oleic acid into respective ester. The results showed that the conversion rate was down with the increase of carbon chain of alcohol, which indicated that high bulk hindrance occurs on the hydroxyl of the alcohol, and the efficient attack of them to the polarized carbonyl of oleic acid is reduced. However, it is not clear how the carbonyl is polarized by Sn 2+ .
Kinetic data of esterification of ethanol and oleic acid catalyzed by SnCl 2 (n ethanol : n oleic acid : n SnCl2 = 120 : 1 : 0.01) were measured at the reflux temperature, and kinetic model was obtained as follows, where C is the concentration of oleic acid, t is the reaction time.
The effectiveness of the catalyst SnCl 2 •2H 2 O has been investigated in a broad range of concentrations, and the results are approximately concomitant with a first order dependence in relation to the catalyst concentration. The effect of the temperature on the initial rate of the esterification was determined, showing that an increase in the reaction temperature was caused a corresponding improvement on the reaction rate, especially at a range of 45-75°C. The value of activation energy for the reaction was determined from the data of the initial rate to be 46.69 kJ·mol -1 .
The kinetics of catalytic esterification of castor oil with lauric acid using SnCl 2 •2H 2 O was studied (Kulkarni & Sawant, 2003 The consumption rate of A (-r A ) is expressed by, where C A , C B , C C is the concentration of lauric acid, hydroxyl group and catalyst (mol/mL), respectively.
The catalyst concentration was constant throughout the reaction; therefore, Equation (11) can be rewritten as The integrated form of Equation (12) using fractional conversion is, where X A , X B is fractional conversion of the lauric acid and the -OH group, respectively; C A0 , C B0 is initial concentration of the lauric acid and hydroxyl group (mol/mL), respectively; and k is the second-order rate constant (mL 2 mol -2 min -1 ).
The esterification reaction of castor oil and lauric acid was carried out using SnCl 2 •2H 2 O as catalyst (0.25,0.5,1.0 and 2.0 % w/w catalyst loadings) at 185℃ and the kinetic data were measured. The results showed that there was a good linearity relationship between ln [(1-X B )/(1-X A )] and t. The plot of k′ versus C C was close to a straight line. For the given values of C B0 , C A0 and C C , the kinetic model of esterification of castor oil and lauric acid was obtained at 185°C, It was investigated the esterification of octanoic acid and n-octyl alcohol utilizing metallic chlorides (KCl, CoCl 2 , MgCl 2 , ZnCl 2 , FeCl 3 etc.) in a stirred tank reactor (Santos, 1996). The results showed that the best efficiency of the formatted ester (n-octyl octanoate) was obtained with ferric chloride, indicating that the higher electronegativity of the metallic ion and the existence of free d orbitals in the transition metals are responsible for the higher yield of ester. Thus, the existence of free d orbitals in Fe 3+ (its configuration: 3d 3 ) gives the possibility of the formation of complexes with OH groups of the reactants, which explains the highest activity found for FeCl 3 •6H 2 O in this way. Generally, in the catalytic systems whose bases were constituted by a transition metal, such as Fe, Co, Mn, Zn, the mechanism of esterification can similarly be ascribed to the above mentioned.
The esterification reaction of octanoic acid and n-octyl alcohol was carried out using CoCl 2 •2H 2 O as catalyst (0, 0.0385, 0.077 mol/l) at 70℃ and the kinetic data were measured. The experimental curves suggest that the kinetics of esterification between octanoic acid and n-octyl alcohol can be described by an irreversible second order power model, considering the catalyst concentration as a constant in the kinetic model proposed. The activation energy is seen to have a value of 53 kcal/mol (Urteaga et al., 1994).
The reaction rate went on a power law model of second order, one for the alcohol and one for the acid, It was also observed that the homogeneous reactions were fitted much better than the pseudo-homogeneous reactions to proposed kinetics model. The reason for it is that an important influence of the physical steps that could have place in the global mechanism of this reaction.

Perchlorate
Perchlorates are of great chemical interest and importance. The high electronegativity, together with the relatively low charge density, results in poor complexing ability of the perchlorate ion. Metal perchlorates can therefore act as powerful Lewis acids, with this character mainly being exploited to activate bidentate compounds (Bartoli et al., 2007). Magnesium perchlorate is one of the most active Lewis acids for esterification. By 2003, Gooβen & Döhring synthesized various esters through a decarboxylative esterification of alkyl and aryl carboxylic acids with dialkyl dicarbonates in the presence of 1 mol% of Mg(ClO 4 ) 2 . However, more in-depth investigations on the mechanism of this reaction are underway. nucleophilicity. Herein, the reaction mechanism was presented to explain the observed reaction characteristics. In fact, the more acidic and less nucleophilic phenols react faster than aliphatic alcohols. Thus, the release of the alcoholic proton should be involved in the ratedetermining step. A reasonable mechanistic hypothesis is depicted as follows (Figure 6), Due to the ability of coordinating with 1,3-dicarbonyl compounds, Mg(ClO 4 ) 2 reacts with diethyl dicarbonate to form complex (II), which can undergo the addition of the alcohol to form intermediate (III). An internal proton shift in (III) can produce intermediate (IV), which can irreversibly decompose to the mixed carbonate (VI) and to the carbonic acid monoester (V). Owing to its high instability, (V) immediately produces EtOH and CO 2 . The irreversibility of the last two steps drives the overall process towards (III). This explanation accounts for the formation of the mixed carbonate (VI) as the major product of the reaction (Jousseaume et al., 2003).
Nevertheless, this is a speculative hypothesis; more studies are in progress to find experimental evidence to elucidate the reaction mechanism.

Organometallic compound
Organometallic compounds containing bonds between carbons and metals provide a source of nucleophilic carbon atoms which can react with electrophilic carbon to form a new carbon-carbon bond (Banach et al., 2001;Finelli et al., 2004). Organometallics find practical uses in catalytic processes. For example, Titanium tetrabutoxide, Butylhydroxyoxostannane are well-known catalysts for the esterification between diacid and diol in polyester industry (Grzesik et al., 2000;Prabhakarn et al., 2011). The catalytic mechanism and kinetic model of related catalysts are reviewed in the following discussions.

Titanate
Titanium-based catalysts have been known for many years and actually are widely used for the production of poly (butylene terephthalate) (PBT), poly (trimethylene terephthalate)

www.intechopen.com
Chemical Kinetics 264 (PTT), and so on. Recently, the developments of titanium-based catalysts make a marvelous progress , such as titanium dioxide based catalyst (C-94) mainly focusing on designing to be stable and having good activity and color.
The mechanism of the mono-esterification between terephthalic acid (TPA) and 1,4-butanediol (BDO) catalyzed by Ti(OBu) 4 was proposed (Tian et al., 2010). As shown in Figure 7, the reaction involves the formation of an adduct between a carbonyl group and Ti atom. Fig. 7. Mechanism of the mono-esterification between BDO and TPA catalyzed by Ti(OBu) 4 Although the reaction system of TPA and BDO is heterogeneous, it can be assumed that the esterification occurs only in the liquid phase. The initial rate method is used to predict the reaction rate. The kinetic model of mono-esterification between TPA and BDO catalyzed by Ti(OBu) 4 in the temperature of 463-483K was investigated (Bhutada & Pangarkar, 1986). The reaction rate r can be described as, where C A , C B is the concentrations of TPA and BDO, respectively, k is the reaction rate constant, n and m is the reaction order with respect to A and B.
The reaction order can be obtained from experiments in which the initial rates are measured at a series of initial reactant concentrations. The results show that the reaction order for TPA and BDO is nearly a constant of 0.7 and 0.9, respectively. Therefore, the rate equation of the esterification is written as follows, where k Ti is the rate constant (L 0.6 ·mol −0.6 ·min −1 ), R 2 is the correlation coefficient.
The reaction of behenic acid with fatty alcohols (decanol, lauryl alcohol, myristyl alcohol and cetyl alcohol) was studied by Tiwari, et al.    (20)) and the catalyzed reactions (Equation (21) (21) Excess of the acids is used to get almost complete conversion of the alcohols, and the unreacted acid was easily removed from the ester as a sodium salt. Assuming that it is first order dependence of the reaction rate on each reactant and the catalyst, the overall rate expression in the integrated form becomes, where X A is fractional conversion of behenic acid; C A0 is initial concentration of behenic acid;  cc u c kk C k, k c is catalyzed reaction rate constant, k uc is uncatalyzed reaction rate constant; C C is catalyst concentration.
The kinetic data show that there is a reasonably good agreement between experimental points and the ones calculated by Equation (22). The values of activation energy obtained for the uncatalyzed (∆E uc ) and catalyzed (∆E c ) reaction of behenic acid with decanol, lauryl alcohol, myristyl alcohol and cetyl alcohol are shown in Table 2.The model Equation (22) is also appropriate for the reaction of erucic acid with cetyl alcohol and oleyl alcohol using TBT as a catalyst

Butylhydroxyoxo-stannane
Stannum-based catalysts such as butylhydroxyoxo-stannane (BuSnOOH) are commonly used for the synthesis of poly (ethylene terephthalate) (PET) , poly (butylene terephthalate) (PBT) and so on. Due to the similarity in coordination number and electronegativity between Ti and Sn, as shown in Figure 8, the mechanism of BuSnOOH catalyst is also proposed by TIAN, et al.
The initial rate method will also be used in this reacting system. The kinetic model of monoesterification reaction between terephthalic acid (TPA) and 1,4-butanediol (BDO) is expressed as follow, where C A and C B is the concentrations of TPA and BDO, respectively; k Sn is the rate constant (L 1.9 ·mol −1.9 ·min −1 ) .

Fluorous metaloxide
In 2004, Otera investigated the transesterification and esterification of various acids and alcohols using 1, 3-disubstituted tetraalkyldistannoxanes as catalyst, and found that the catalytic mechanism different from that proposed for tetraalkyldistannoxanes. In the latter case, the initial step is substitution of the bridging X group to give an alkoxydistannoxane intermediate, (YR 2 SnOSnR 2 OR) 2 , which works as an alkoxyl donor. In contrast, for 1, 3disubstituted tetraalkyldistannoxanes, no substitution takes place at the bridging position b y a n a l k o x y g r o u p . A s d e s c r i b e d a b o v e , t h e b r i d g i n g c h l o r i n e i n 1 , 3 -d i s u b s t i t u t e d tetraalkyldistannoxanes is never substituted by isothiocyanate ion, but experiences strong association.
A related mechanism of the esterification is proposed that both alcohol and acid coordinate on the terminal tin atom (Yoshida et al., 2006), on which the interchange between the hydroxyl groups of carboxylic acids and the alkoxy groups of alcohols takes place. Meanwhile, since water, one of the products, is less fluorophilic, the esterification of carboxylic acids with alcohols should proceed efficiently.
Fluorous distannoxanes exhibit unusually high preference for fluorous solvents over common organic solvents thanks to coverage of the molecular surface with fluoroalkyl groups. While the fluorous biphase technology had been high-lighted mainly in terms of facile separation of products and fluorous catalysts, it has now been revealed that the equilibrium is also controllable under fluorous biphase conditions. Thus, the esterification has been completely driven in the desired direction by use of a 1:1 ratio of starting materials without recourse to any dehydration technique.

Solid acids
The research for solid acids has become active since the early 1970s. In 1971, Isao, et al. investigated the esterification of ethanol with acetic acid on silica-alumina, and a simple kinetic model based on a Langmuir-Hinshelwood mechanism was proposed. Since then, the esterification catalyzed by solid acids is widely studied and largely reported (Jiang et al., 2008;Jothiramalingam & Wang, 2009;Li et al., 2010). Herein, the esterification mechanisms and kinetic models according to the different types of solid acids are reviewed, such as sulfate-supported metal oxides (SO 4 2− /M x O y ), TPA/SnO 2 .

TPA/SnO 2
The esterification of palmitic acid with methanol using 12-Tungstophosphoric acid (TPA)/SnO 2 was investigated (Srilatha et al., 2011). The reaction is shown as follows: As shown in Figure 9, the mechanism of this esterification is proposed,  The catalyst initiates the esterification reaction by donating a proton to palmitic acid molecule. The palmitic acid is then subjected to nucleophilic attack by the hydroxyl group of methanol, and the reaction continues with water elimination.
As the large excess of methanol, it could be safely assumed to be a first-order pseudohomogeneous reaction. The esterification reaction was carried out using 15 wt% TPA/SnO 2 as catalyst w A : w cat =5:1 in the temperature range of 45-65℃and the kinetic data were measured.
The reaction rate can be described as follows: where X A is the conversation of palmitic acid, t is the reaction time, k is the reaction rate constant.

Chemical Kinetics 268
The plots of-ln (1-X A ) versus t at different temperatures showed that there was a good linearity relationship. Linear fitting to lnk -1/T curve was carried out and the reaction rate constant can be described as follows: The activation energy was obtained as 36.33 kJ mol −1 .

SO 4 2− /M x O y
The solid superacid (SO 4 2− /M x O y ) is a new type of catalyst used in esterification (Jiang et al., 2004). M x O y are usually some transition metal oxides such as ZrO 2 and TiO 2 .
In 2010, Rattanaphra, et al. investigated the esterification of myristic acid with methanol catalyzed by sulfated zirconia. The mechanism of esterification is probably following Langmuir-Hinshelwood model (Arata, 2009;Reddy & Patil, 2009). It is possible that methanol and myristic acid are preferentially adsorbed on the Bronsted acid sites of sulfated zirconia during esterification. The hydroxyl group of methanol is protonated by Bronsted acid on the catalyst surface while the protonation of myristic acid on an adjacent site leads to the carbocation. Deprotonation of methanol oxygen produces the nucleophile, which attacks the carbocation to generate a tetrahedral intermediate. As shown in Figure 10 The stirring rate is sufficient to overcome the diffusion limitation of reactive species. Therefore, the performance of pseudo-homogeneous model can be considered as satisfactory to correlate the kinetic data for the esterification. The rate equation can be written as follows: where C A , C B , C E and C W is the concentration of myristic acid, methanol, myristic acid methyl ester and water respectively, k 1 is the forward rate constant and k -1 is the backward rate constant.
The equation can be rearranged to be: where X A is the conversion of myristic acid，C A0 is the initial concentration of myristic acid, K e is the equilibrium constant, M is the concentration ratio of methanol to myristic acid (M=C B0 /C A0 ).

   
The linear coefficient is 0.996, and the activation energy is 22.51 kJ /mol.
In 2004, Jiang, et al. studied the esterification of n-pentanol with benzoic acid using Alpillared clay (PILC) supported SO 4 2− /TiO 2 superacid catalyst (Fang et al., 2010), and find that it is known that the generation of super acid sites in the system of SO 4 2− /M x O y solid superacid is necessarily promoted by the sulfur of the metal oxides, the more acid sites formed, the higher catalytic activity exhibited. Therefore, Al-PILC carrier can effectively enhance the catalytic activity.
In 2010, Zubir, et al. investigated the kinetic behavior of the heterogeneous catalyzed esterification of oleic acid with ethanol using tungstated zirconia as a catalyst in the temperature range of 303.15-323.15K (Otera, 1993 The kinetics of esterification reaction can be expressed using a pseudo-homogeneous second-order equilibrium model in the absence of any intraparticle diffusional limitation as follows: www.intechopen.com

270
where C A ,C B , C C , C D is the concentration of oleic acid, ethanol, ethyl oleate and water, respectively, k is the kinetic constants for the forward reaction, K e is the equilibrium constant.
The value of activation energy is 51.9 kJ·mol -1 . The goodness-of-fit of the experimental data to the proposed model is assessed by comparing the experimental reaction rate with the theoretical prediction, and the experimental data are reproduced with errors not greater than 10%.

Ion-exchange resin
Ion-exchange resins, especially the cation-exchange resins such as Dowex, Amberlyst series are manufactured mainly by sulfonation of ethylbenzene first, followed by a cross-link with divinylbenzene (Liu & Tan, 2001;Alexandratos, 2008;Tesser et al., 2010). Because of their selective adsorption of reactants, surface acid site features, and swelling nature, these resins not only catalyze the esterification reaction but also affect the equilibrium conversion. They also show excellent performance such as reusable, mechanical separation, continuous operation as a heterogeneous catalyst in esterification (Yang et al., 2007;JagadeeshBabu et al., 2011;Ju et al., 2011;.
By 1965, Bochner, et al. evaluated the performance of Dowex 50WX-8 as a catalyst for the esterification of salicylic acid with methanol. Through analyzing the experimental results, the Langmuir-isotherm applied to this system to describe the reaction rate. Since then, various mechanisms and kinetics of esterification between acids and alcohols catalyzed by cation-exchange resins were proposed. Here is an overview of the reaction mechanisms and kinetic models according to the type of cation-exchange resins.

Dowex
Dowex, a common cation-exchange resin, is widely used as catalysts in esterification (Vahteristo et al., 2009). In 2005, Kulawska, et al. investigated the esterification of maleic anhydride with octyl, decyl or dodecyl alcohol over Dowex 50Wx8-100, and the kinetic data were measured in the temperature range of 403-433 K (C 0ALC /C 0MA =5:1). Based on the data, a first order reaction was found-first order with respect to acid and zero order with respect to alcohol. The reaction rate can be described as follows, where c M is the concentration of maleic monoester, t is reaction time.
The values of activation energy are 66.0 (±0.65) kJ/mol, 58.6 (±0.4) kJ/mol, and 66.1(±0.4) kJ/mol for dioctyl, didencyl, and didodecyl maleate formation, respectively. However, the effects of the catalyst (the particle size and the concentration) and the mass transfer resistance on the reaction rate were not taken into account.

Recent Developments on the Mechanism and Kinetics of Esterification Reaction Promoted by Various Catalysts 271
The esterification kinetics of acetic acid with iso-butanol catalyzed by Dowex 50 Wx2 was studied (Izci & Bodur, 2007). It is considered to be reversible reactions and can be described by pseudo-homogeneous (PH) model. The general reaction rate expression can be written as follows, where subscripts A, B, E and W is acid, alcohol, ester and water, respectively, k 1 is forward reaction rate constant (L mol −1 min −1 ), K is the equilibrium constant of the reaction.
It is found that this bimolecular type is second order reaction. The values of k 1 for different temperatures are 0.24×10 3 (318K), 1.00×10 3 (333K), 3.07×10 3 (348K). The activation energy was found to be 1.745 kJ·mol −1 in the presence of Dowex 50 Wx2. However, the PH model does not take into account the resin swelling ratio, adsorption of the components and the non-ideal thermodynamic behavior of reactants and products (Ali &Merchant, 2006;Patel & Saha, 2007;Jeřábek et al., 2010;Erdem & Kara, 2011).
The kinetics of Dowex 50 Wx8-catalyzed esterification was studied between acetic acid and benzyl alcohol (Ali & Merchant, 2009). The swelling ratio of Dowex 50 Wx8 in different solvent was measured, and the results show that it decreases in the order of water, benzyl alcohol, acetic acid, benzyl acetate. Water appears to be preferentially adsorbed by the catalyst from a binary solution of acetic acid and water, and hinders the approach of butanol to the protonated acid. Therefore, the water exerts an adverse effect on the esterification rate.
The initial reaction rates of esterification between acetic acid and benzyl alcohol were measured at various conditions, and the Eley-Rideal (ER) model was used to correlate the data and showed a high degree of fit, indicating that the surface reaction between adsorbed alcohol and acid in the bulk is the rate-limiting step during the initial stage of the reaction.
where M cat is the mass of the catalyst, g; a acid , a alc , a ester , a water is the activity of acid, alcohol, ester, water in the liquid phase, respectively; k f is forward reaction rate constant for the esterification, mol g -1 s -1 . K a is activity reaction equilibrium constant. K acid , K alc , K water is adsorption equilibrium constant for the alcohol present in the system, respectively.
The activation energy form the above relationship was found to be 55.4kJ/mol.
However, as the reaction proceeds, the acid adsorption term might have to be introduced and in such an event the reaction kinetics would be represented by a dual site mechanistic, such as Langmuir-Hinshelwood (LH) model.
In 2011, Ju, et al. measured the kinetic data for the esterification of butyric acid with nbutanol over Dowex 50Wx8-400, and correlated with various types of kinetic models. Strong resin water affinity was taken into account, and the non-ideality of the system was considered by applying the no-random two liquid (NRTL) model. The comparisons of the conversion of butyric acid with reaction time between experimental data and ones predicted by these kinetic models reveal that the ER model and LH model are the most reasonable fit for describing the mechanism, with the total average error 12.52%. Surface reaction is the rate determining step, and the affinity between resin and water is found to be not strong. Therefore, the most possible esterification reaction mechanisms can be proposed as following lists: 1. Single site mechanism: the adsorbed butyric acid onto the catalysts reacts with nonadsorbed n-butanol in the bulk. 2. Single site mechanism: the adsorbed n-butanol onto the catalyst reacts with nonadsorbed butyric acid in the bulk. 3. Dual site mechanism: both the reactants adsorbed on the catalyst surface and react there.
In 2007, Ali, et al. studied the esterification of 1-propanol with propionic acid catalyzed by Dowex 50Wx8-400. The experiments were carried out over a temperature range of 303.15-333.15K, and the reaction mechanism for the esterification was proposed (Lilja et al., 2002). The reaction is initiated by the transfer of a proton from the catalyst to the carboxylic acid，and the carbonium ion is formed during the reaction. The ion is accessible for a nucleophilic attack by the hydroxyl group from the alcohol. After that, a molecule of water is lost from the ion. Finally, the catalyst is recovered by the transfer of proton from the ion to the catalyst surface. This mechanism is represented by the following scheme ( Figure 11): Fig. 11. Mechanism for the esterification catalyzed by Dowex 50Wx8-400 The donation of a proton is commonly assumed to be a fast step, and the nucleophilic substitution is usually assumed to be the rate determined step.
Based on the above reaction mechanism, external and internal diffusion can be negligible, and the E-R model (Equation (39)) is applied. The activity coefficients of components can be predicted using UNIFAC model. Meanwhile, by taking into account the strong water affinity for Dowex 50Wx8-400, a correction term (α) was added to the activity term for water in rate expression. The experimental results for the esterification between 1-propanol and propionic acid were modeled according to the above kinetic model. The total average error between the predicted and experimental mole fractions of acid is 1.65%. The activation energy for the forward reaction was estimated to be 67.3 kJ /mol.
The effect of different alcohols on the conversion of propionic acid was investigated using methanol, ethanol, 1-propanol and 1-butanol. The results reveal that the conversion of The acetic acid initiates swelling of the resin, which results in easy accessible acid groups for the reaction and free mobility of all the components. Therefore, a pseudo-homogeneous (PH) kinetic model is applicable. For the reaction, there is, where k 1 , k 2 , k −1 , k −2 is the forward and backward rate constants, respectively, (mol g -1 s -1 ), a i is the activity of each component, and it is calculated using the UNQUAC model, m cat is the mass of catalyst.
The kinetic behavior of esterification of lactic acid with isopropanol over Amberlyst 15 was investigated Kirbaslar et al., 2001), at different temperatures from 323K to 353K. The ER model was used to describe the reaction mechanism that takes place between adsorbed molecules of isopropanol and the molecules of lactic acid in bulk. The adsorption of ester is reported to be negligible. Hence the rate equation can be described as follows, where k 1 is Rate constant for esterification reaction, M is molar ratio of isopropanol to lactic acid, K e is equilibrium rate constant.
The value of activation energy is found to be 22l J/mol, and at each different temperature, the rate constants is 1. The heterogeneous esterification of propionic acid with n-butanol over Amberlyst 35 for the synthesis of n-butyl propionate was investigated (Lee et al., 2002). The esterification is shown as follows,   The kinetics data of the liquid-solid catalytic esterification are correlated with various kinetic models, over wide ranges of temperature and feed composition. The activity coefficients calculated using the NTRL model are utilized to represent the non-ideality behavior of the species in the liquid solutions. Meanwhile, the effects of film diffusion and pore diffusion appear to be negligible at the experimental conditions. The results reveal that the Langmuir-Hinshelwood (LH) model yielded the best representation for the kinetic behavior of the liquid-solid catalytic esterification: The activation energy of reaction is about 6.81kJmol −1 . This model is capable of representing the kinetic behavior of the liquid-solid catalytic esterification at temperatures from 353.15 to 373.15 K over entire range of the experimental feed compositions.
The kinetics of esterification between acrylic acid and propylene glycol in the presence of Amberlyst 15 was investigated (Altıokka & Ödeş, 2009). Taking into account the general esterification reaction as well as polymerization of acrylic acid and products, the overall reaction mechanism is proposed to be: The reaction rate constants k 1 , k 2 , k 3 and k 4 , in Equation (54) are determined, and shown in Table 4, where k 1 , k 2 is forward and backward reaction rate constants in Equation (51), and k 3 , k 4 is forward reaction rate constants in Equation (52) and Equation (53), respectively,(L min −1 mol −1 ), T is absolute temperature in K.
Together with the experimental data, the concentration-time curves based on the model were obtained under given reaction conditions; there is a reasonably good agreement between calculated curves and experimental points.

Conclusions
This chapter discussed esterification mechanisms, and evaluated the kinetics objectively and quantitatively, which provided a most effective way to select catalyst and design reactor for different esterification systems. It is discovered that some new catalysts (such as lipases, room temperature ionic liquids) have being used in esterification; nevertheless, there are few research on the case. Herein, it is worthy to be investigated deeply.


extfieldsofstudy:
Medicine
Computer Science

text:
An integrated open framework for thermodynamics of reactions that combines accuracy and coverage

Motivation: The laws of thermodynamics describe a direct, quantitative relationship between metabolite concentrations and reaction directionality. Despite great efforts, thermodynamic data suffer from limited coverage, scattered accessibility and non-standard annotations. We present a framework for unifying thermodynamic data from multiple sources and demonstrate two new techniques for extrapolating the Gibbs energies of unmeasured reactions and conditions. Results: Both methods account for changes in cellular conditions (pH, ionic strength, etc.) by using linear regression over the ΔG○ of pseudoisomers and reactions. The Pseudoisomeric Reactant Contribution method systematically infers compound formation energies using measured K′ and pKa data. The Pseudoisomeric Group Contribution method extends the group contribution method and achieves a high coverage of unmeasured reactions. We define a continuous index that predicts the reversibility of a reaction under a given physiological concentration range. In the characteristic physiological range 3μM–3mM, we find that roughly half of the reactions in Escherichia coli's metabolism are reversible. These new tools can increase the accuracy of thermodynamic-based models, especially in non-standard pH and ionic strengths. The reversibility index can help modelers decide which reactions are reversible in physiological conditions. Availability: Freely available on the web at: http://equilibrator.weizmann.ac.il. Website implemented in Python, MySQL, Apache and Django, with all major browsers supported. The framework is open-source (code.google.com/p/milo-lab), implemented in pure Python and tested mainly on Linux. Contact: ron.milo@weizmann.ac.il Supplementary Information: Supplementary data are available at Bioinformatics online.

aims to better understand the complex behaviour of living organisms as well as their manipulation for human needs (Atsumi et al., 2008;Bar-Even et al., 2010;Sinha et al., 2010;Steen et al., 2010;Yim et al., 2011;Zhang et al., 2010). The thermodynamics of biochemical reactions (Alberty, 2003) is of special interest in the analysis and design of metabolic pathways.
The change in the Gibbs free energy ( G) characterizes the thermodynamic balance of biochemical reactions and dictates the direction of net flux (the difference between forward and backward fluxes) in a reaction. It is thus useful for the study of a single enzymatic reaction, for analyzing entire metabolic pathways (Vojinović and von Stockar, 2009), and for the large-scale modeling of whole-cell metabolic networks (Henry et al., 2007).
A reaction at equilibrium carries no net flux. At equilibrium in a specific pH, the apparent reaction quotient Q -the ratio of product to substrate concentrations-is termed the apparent equilibrium constant and denoted by K . In ideal dilute solutions, the 'transformed Gibbs energy of reaction' is a function of the apparent reaction quotient: r G =−RTln K /Q . The 'standard' transformed Gibbs energy of reaction ( r G • ) is the value of r G at standard conditions, i.e. when all compound concentrations are 1 M. Therefore, r G • =−RTlnK ( Fig. 1-equation I).
r G determines the direction of net flux in a reaction. A negative r G would correspond to a positive (forward) net flux and vice versa. Cell physiology imposes constraints on metabolite concentrations and consequently on Q . Reactions for which r G • < 0 for any physiological Q can only carry a forward flux and are thus called irreversible reactions. This classification of reactions is especially important in constraint-based modeling that covers whole-cell metabolic networks and depends on the knowledge of reaction directionality for predicting flux distributions, growth rates and other large-scale metabolic phenotypes (Beg et al., 2007;Burgard et al., 2003;Oberhardt et al., 2009). Directionality, annotations typically rely on phenomenological data and arbitrary definitions of reversibility. Recent advances in the field allow incorporating thermodynamic data directly into the model by adding explicit constraints that connect r G, concentrations and reaction directionality (Fleming et al., 2009;Henry et al., 2007).
The NIST database for Thermodynamics of Enzyme-Catalyzed Reactions (NIST-TECR) is the most comprehensive collection of empirical thermodynamic data (Goldberg et al., 2004(Goldberg et al., , 2007. About 400 reactions have measured equilibrium constants and were successfully mapped to KEGG identifiers (Kanehisa et al., 2008). Still, this represents only about 6% of the ∼5500 biochemically relevant reactions in the KEGG database. Moreover, since equilibrium constants in the NIST-TECR database are measured in a variety of pH and ionic strength levels, it is difficult to extrapolate the equilibrium constant for the conditions prevailing in an organism or system of interest. When r G • are not covered by the NIST-TECR database, several collections of standard Gibbs free energy of formation ( f G • ) can be consulted and used as in Figure 1equation II. The largest collection, given by Alberty (2003), contains f G • values for ∼200 compounds. Unfortunately, many reactions were not measured in the desired conditions (pH, ionic strength, etc.) or have never been experimentally measured at all. If the formation energy of even one of the reactants is unknown, r G • of a reaction cannot be derived. In order to bridge the gap between the known formation energies and the unknown ones, a method based on the group contribution assumption for biochemical compounds in aqueous solutions was described by Mavrovouniotis (1990Mavrovouniotis ( , 1991 and later greatly improved in terms of coverage and accuracy by Jankowski et al. (2008). The method is based on the simplifying assumption that each functional group has a characteristic contribution to the overall formation energy of a compound and that these group contributions are independent of each other. Therefore, f G • is estimated by summing all contributions from the different groups composing a compound. The contribution of each such group ( gr G • ) is estimated through linear regression, which uses the known f G • and r G • and the compounds' partition into groups (Mavrovouniotis, 1990).
The group contribution method is limited in its accuracy. The group independence assumption can be overly simplistic, especially for large compounds or in conjugated systems. In addition, the definitions of the groups rely heavily on heuristics and chemical intuition. Recent improvements of the method included a better choice of groups and the introduction of group interaction corrections (Jankowski et al., 2008). Recently, a promising new approach based on whole-reaction similarities (Rother et al., 2010) has been shown to be more accurate, but currently does not provide as wide a coverage as group contribution methods.
Previous group contribution studies did not consider the effect of pH on the compounds' protonation levels (Alberty, 1998;Jankowski et al., 2008;Mavrovouniotis, 1990Mavrovouniotis, , 1991 and assumed all measurements were taken at standard aqueous conditions (e.g. pH 7 and ionic strength of 0.25 M). Each compound is actually an Fig. 1. An overview of the relationships between the layers of thermodynamic data. Reaction equilibrium constants can be converted to r G • which is calculated as the stoichiometric sum of reactant formation energies (blue shading, equations I and II). Each reactant exists as several pseudoisomers distributed according to the Boltzmann distribution, and its f G • can be found using the Legendre transform (yellow shading, equations III and IV). Using the group contribution assumption, each pseudoisomer can be constructed from its group components (green shading, equations V and VI) ensemble of 'pseudoisomers' differing in their protonation levels (e.g. ATP 4− , HATP 3− and H 2 ATP 2− ) and its f G • can be found using the Legendre transform ( Fig. 1-yellow shading, equations III and IV). As pointed out by Alberty (2003), ignoring the change in the abundance of the different pseudoisomers at changing pH levels can result in errors. This is especially true for the biochemically ubiquitous phosphate groups, which have pK a values in the range of 6-8 (Lide, 2009;Robitaille et al., 1991), and hence their pseudoisomer distribution changes considerably even close to pH 7.
The concentrations of the individual pseudoisomers are determined by the Boltzmann distribution (Alberty, 2003). f G • (i) of a compound i is, therefore, a function of f G • ( j) of the compound's pseudoisomers ( j) (Fig. 1-equation III). Notably, f G • ( j) of the pseudoisomers and their distribution is modulated by the pH (Fig. 1-equation IV). For example, at a high pH, the unprotonated form of an acidic compound will have lower f G • and thus be more abundant (Alberty, 2003). In order to use f G • to calculate r G • , the different pseudoisomer forms assumed by each compound should be considered.

APPROACH
In this study, we present methods to accurately estimate r G • using two major approaches (Fig. 1). The first approach, named Pseudoisomeric Reactant Contribution (PRC), recovers pseudoisomer formation energies by applying linear regression to the entire set of reactions available in the literature (as recorded in NIST-TECR). Most measurements in the NIST-TECR database are of 'apparent' equilibrium constants, which is the equilibrium apparent reaction quotient of the total compound concentrationsi.e. the sum of all protonation states. Since the pH affects the distribution of species non-linearly, the value of r G • cannot be expressed as a linear sum of the reactants' formation energies. PRC applies the inverse Legendre transform to linearize the system. The inverse Legendre transform replaces each ensemble of pseudoisomers with a single representative and changes the observed value of K accordingly (Alberty, 2002) using the known dissociation constants of the reactants. It is thus possible to find the least-squares solution for the f G • of the pseudoisomers using linear regression (see Supplementary Material for details). Using PRC, we can retrieve f G • values that were previously unknown (Alberty, 2003) and enable the calculation of r G • for more reactions.
Compounds that were not measured previously and thus do not appear in the NIST-TECR database require a different approach. We have developed an augmented group contribution method that we call Pseudoisomeric Group Contribution (PGC). Unlike previous approaches, we decompose pseudoisomers, not compounds, into functional groups ( Fig. 1-green shading). We then estimate f G • of the pseudoisomers by summing over the contributions of their groups (equation V) and calculate the f G • of compounds (equations III and IV). This method provides higher accuracy and can correctly adjust the f G • to the different aqueous conditions described by pH and ionic strength (with pH usually being the most significant). The combination of the two approaches enables a better estimation of f G • for a large variety of compounds in a wide range of aqueous conditions.

A PRC method systematically derives f G •
The task of estimating the formation energies of compounds given measured reaction equilibrium constants is not straightforward. The difficulty stems from the non-linearity of r G • as a function of f G • and pH ( Fig. 1equations II-IV). However, if r G • has been measured in a specific pH and all but one of the reactants' f G • values is known, it is straightforward to infer the missing f G • .
The extensive list of compound formation energies provided in (Alberty, 2003) is the product of a meticulous reconstruction such as described above and based on the data provided by many measurements of K . As each new f G • added to this database relies on the previously gathered values, measurement errors for a particular compound are carried on to future calculations.
Alternatively, compound acid dissociation constants (pK a s) can be used to convert the set of equations into a linear system that can be solved computationally. The idea, known as the 'Inverse Legendre Transform' (see Supplementary Material), is based on the fact that the difference between the r G • and r G • of any reaction is a function of the pK a s alone and does not depend on the absolute formation energies of the reactants. The resulting linear system can then be solved using multiple linear regression. This application of the inverse Legendre transform to linearize a reaction system was introduced by Alberty (1991) and Alberty and Goldberg (1992), but to our knowledge was implemented only for a small set of reaction measurements. Here, we perform a global analysis using all the available data in NIST-TECR to achieve the best possible least-squares estimation. We refer to this method as the Pseudoisomeric Reactant Contribution (PRC) method.
Using PRC, we were able to obtain values for 407 formation energies using only the 367 reactions in the NIST-TECR dataset and the pK a values of the participating compounds (described later in 'Acid dissociation constants'). A detailed analysis of these values in terms of prediction power and accuracy is given in Section 4. Due to linear dependencies between some of the reactions, there are 112 dimensions in the null space of the stoichiometric matrix. For instance, if two compounds always appear together (such as NAD ox and NAD red ), the difference between the f G • of the pair can be inferred, but the absolute formation energy remains unknown. Similarly, element conservation rules contribute one dimension to the null space for every element which appears in the database, namely C, N, O, S and P. Commonly, as is the case in the current study, the ambiguity in the values of the formation energies is solved by defining some compounds as having f G • = 0. Alberty's table of ∼200 formation energies (Alberty, 2003) contains 18 such reference points.

A PGC method covers more reactions and conditions
In order to improve the estimations provided by the group contribution method (Mavrovouniotis, 1990) for highly pH-dependent compounds, we introduce a method that incorporates Alberty's transformed formation energies (Alberty, 2003) into the same framework used by Jankowski et al. (2008). Previous group contribution implementations used apparent equilibrium constants (K ) and formation energies ( f G • ) at pH ∼7, where each compound is actually an ensemble of protonation species-or pseudoisomers ( Fig. 1-equation III). For instance, the total concentration of ATP is divided between the pseudoisomers H 2 ATP 2− , HATP 3− and ATP 4− according to the Boltzmann distribution. The formation energy of the ensemble is called the standard transformed Gibbs energy as defined by the International Union of Biochemistry and Molecular Biology (IUBMB) (Alberty et al., 2011). A shift in pH will change the relative abundance of the different pseudoisomers and affect f G • and K nonlinearly. Therefore, in order to normalize the effect of pH across different measurements, we use only formation energies of single pseudoisomers as input for the linear regression step in our group contribution framework. For example, instead of the 'transformed' formation energy of ATP (-2292.5 kJ/mol), we use the formation energy of HATP 4− , which is -2768.1 kJ/mol (Alberty, 2003). This change requires knowing the standard formation energy of each pseudoisomer separately and the exact distribution of protons and charges in its molecular structure. In the case of reactions, we use the inverse Legendre transform-similarly to the PRC method.
In addition, groups that previously had only one form corresponding to the protonation level most abundant at pH 7 can now have multiple forms representing the different protonations at a wide range of pHs. The full list of groups is given in Supplementary Table S4. For example, the terminal phosphate group −PO 2− 3 which had only one instantiation in previous implementations (Jankowski et al., 2008) ( gr G • =−254 kcal/mol =−1063 kJ/mol), has two versions in PGC, namely −PO 2− 3 and −HPO − 3 , each with its own contribution ( gr G • =−1024.2 and −1073.9 kJ/mol, respectively). The terminal phosphate of ATP's major pseudoisomer at pH 7 is semi-protonated and corresponds to the −HPO − 3 group. However, when this group appears in the major pseudoisomer of other compounds such as d-glucose-6-phosphate, it is fully deprotonated (−PO 2− 3 ). Thus, PGC adjusts the contribution of the two pseudoisomeric groups accordingly. Similarly to previous implementations (Jankowski et al., 2008), the algorithm for determining the group contributions is a least-squares linear regression, except that for the same number of compounds, there are more available values of f G • and more groups due to the use of pseudoisomers.
We use the NIST-TECR database (Supplementary Table S1), together with a list of formation energies (Alberty, 2003;Dolfing and Janssen, 1994;Thauer et al., 1977) (Supplementary Table S2) and dissociation constants (Lide, 2009) (Supplementary Table S3). As part of this study, we have manually annotated each of the species in the list and determined the protonation level of each of their groups. The task was not trivial for compounds with more than one protonation site, since the order of the pK a values of the groups determines which will deprotonate first as the pH rises, and these data are not well organized in the literature. For example, the different protonation states of nucleic acids are a particular challenge (Christensen et al., 1970a, b). The details and results of this analysis and the contribution of the pseudoisomeric groups are given in Supplementary Table S4.
Using the PRC method to infer compound formation energies from NIST-TECR, it is possible to obtain predictions for a few hundreds of reactions in the KEGG database (roughly 5%)-not much more than the number of reactions in NIST-TECR itself. However, the PGC method achieves a much higher coverage of ∼80% of KEGG reactions, as shown in Table 1. This coverage is similar to previous implementations (Jankowski et al., 2008). Note that although there are >8000 reactions listed in KEGG, we limit our analysis to the subset of reactions whose reactants all have explicit chemical formulas (unlike wildcard formulas or generalized names such as 'donor' and 'acceptor'). Furthermore, we discard reactions that are not chemically or redox balanced.

Acid dissociation constants
Calculator Plugins were used for structure property prediction and calculation, provided by Marvin 5.5.1, 2011 from ChemAxon (http://www.chemaxon.com). The molecular description of each compound in the KEGG database was converted to a SMILES string, and given as input to the ChemAxon command-line binary which calculated the pK a values and the charge distribution of the major pseudoisomer at pH 7. All results are listed in the Supplementary Dataset.

The NIST-TECR database benchmark
In order to evaluate the level of error for each method, we compare the measured r G • and the estimated one (adjusted to the same conditions as the measurement). If the chosen estimation method does not cover this  Jankowski et al. (2008), the reported result is 1.90 kcal/mol (8.0 kJ/mol). In the current article, more reactions from NIST-TECR were used for calculating RMSE since we did not filter observations according to their pH. In addition, we did not include any formation energy data in this evaluation since they are not purely empirical-many are derived from the same data which is already in NIST-TECR, and some values were derived using a group contribution approach and thus cannot be used to evaluate its precision.
reaction, its observed equilibrium constants are not included in this analysis. In order to minimize the temperature related biases, only measurements at temperatures in the range of 298-314 K (i.e. 25-40 • C) were used. The root mean squared error (RMSE) of each method is calculated by giving each distinct reaction an equal weight (regardless of how many times it has been measured). The result for a leave-one-out cross-validation test for the PGC method was 11.2 kJ/mol (RMSE). All formation energy and reactions examples were included in this test. If an example could not be omitted because it was linearly independent of the other examples in the training set, its error was not included in the test. If several examples (reactions or compounds) had exactly the same group decomposition they were all treated as a single example. We used our framework to perform a cross-validation test for the standard (i.e. non-pseudoisomeric) group contribution method, which resulted in an RMSE of 16.0 kJ/mol (see Supplementary Material). It was not possible to directly compare our cross-validation results to the other methods mentioned as their implementations are not publicly available.

The reversibility index
We introduce a quantitative measure for reversibility that takes the mean concentration of metabolites and the number of substrates and products into consideration. Further details regarding the application of the index to genome-scale models and motivation for its use are given in Section 4.
We assume that all reactant concentrations lie within a range located symmetrically (in log-scale) around some characteristic physiological concentration C-here we use 100 μM (Bennett et al., 2009). The variable describing this range is denoted γ : all substrates are assigned a concentration of (1/ √ γ )C and all products have a concentration of √ γ C. Therefore, the log-scale width of the range is log( √ γ C)−log((1/ √ γ )C) = log(γ ). Thus, a value of γ = 100 corresponds to the range 10 μM-1 mM and γ = 1000 corresponds to about 3 μM-3 mM.
For a reaction with N S substrates and N P products, the apparent reaction quotient Q (γ ) = (γ 1/2 C) N P ·(γ −1/2 C) −N S = γ N P /2 γ N S /2 C N P −N S ≡ γ N /2 Q , where we define N = N P +N S as the total number of reactants and Q = C N P −N S as the default reaction quotient. We define the reversibility index asγ = K /Q 2/N , which is the required concentration range for reversing the reaction, i.e. Q(γ ) = K . The furtherγ is from 1 the more irreversible the reaction. The reversibility index of the fructose-bisphosphate aldolase reaction, for example, is 1.04, making it clearly reversible since a change of only 4% in concentrations is required to reverse its direction.

Fig. 2.
A comparison between the apparent equilibrium constant K and the reversibility indexγ . Each shaded square represents a different regime where a reaction is considered reversible or irreversible according to K andγ , and the number indicates the percent of reactions in the metabolism of E.coli (in the iAF1260 model) which occupy that regime. White areas represent regimes where both classifications agree, and pink areas are regimes where they disagree. Colored lines show the relationship between K andγ for specific reactions stoichiometries and the percent of reactions in the model with that stoichiometry is given in parentheses. The inset on the bottom right shows the fraction of reversible reactions (Y -axis) as a function of the allowed concentration range-determined by the upper bound on the value of γ (X -axis). This is equivalent to the cumulative distribution function (CDF) of the reversibility index for all reactions in the model. Note that for the CDF, the direction of each reaction is defined such thatγ will be larger than 1 It should be noted that for reactions with one substrate and one product, Q = 1 and N = 2, and thereforeγ = K . However, since 88% of reactions are not 1 : 1,γ is usually very different from the value of K , like in the case of fructose-bisphosphate aldolase whereγ = 1.04 while K ∼ 10 −4 . Evaluating both K andγ for all the reactions in the E.coli iAF1260 model (Feist et al., 2007) reveals that K plays only a partial role in determining the reversibility index (Fig. 2). In the inset, we show the fraction of reversible reactions as a function of allowed physiological concentration range reflected byγ .
In some cases, it is useful to assign a prescribed concentration for a reactant, rather than allowing it to lie within the range defined by γ . For instance, the concentration of water is fixed and in some organisms, the concentration of P i is kept almost constant. In such cases, the required concentration can be used to calculate Q instead of the default C and the fixed reactant is not included in the value of N .

Source code
The code for the implementation of the PRC and PGC methods is free and open source and can be found at: http://milo-lab.googlecode.com/svn/branches/bioinfo-2012/ Our software is written completely in Python, depending only on free software such as Open Babel (openbabel.org) and SciPy (scipy.org).

Fig. 3.
The response of PPDK to pH, as recorded in the NIST-TECR database. The ionic strength for all estimation methods was set to the default 0.1 M. By taking only the most abundant psuedoisomers at pH 7 and not updating to pK a values, the charge through this reaction is +1, i.e. one proton is added to the products. Therefore, the predicted response to pH will have a positive slope of RTln(10) (blue line). The intercept of the line was determined so that the value at pH 7 would match the predictions of Jankowski et al., 2008. PGC and PRC use pK a values derived by ChemAxon (see Section 3) to calculate the response, which results in a better fit to the measured data (orange and green lines). Note that the lines for PGC and PRC are not exactly straight due to the combined effect of several pK a values. The measured data (dark green) which is listed in NIST-TECR, was taken from (Reeves and Menzies, 1968)

Adjustments for pH using pK a values can increase accuracy
Many organic compounds are weak acids and bases and, as such, assume multiple protonation levels in typical physiological conditions. Formation energies and, as a result, reaction energies are a function of the distribution of reactant psuedoisomers ( Fig. 1equations II and IV). Since the distribution of protonation levels is a function of the prevailing pH and ionic strength, f G • and r G • vary with pH and ionic strength as well. For example, the key reaction in gluconeogenesis and C4 plant photosynthesis, pyruvate + ATP + P i PEP + AMP + PP i , is catalyzed by PPDK (pyruvate-phosphate dikinase, EC number 2.7.9.1) and was measured at various pH levels ranging from 6.5 to 8.4 (Fig. 3). Many of the reactants have a pK a in this pH range, and so r G • of the PPDK reaction changes significantly with pH. Intracellular pH values are typically between 6 and 7.5 for most organisms (Vojinović and von Stockar, 2009), but the range can be much wider, e.g. in yeast (5.5-7.5 (Imai and Ohno, 1995;Ryan and Ryan, 1972)) or bacteria (4.8-9.0 (Breeuwer et al., 1996)). As shown in Figure 3, the PGC and PRC methods can accurately predict the pH dependence of the reaction. Methods that do not account for pseudoisomers can result in a difference of >20 kJ/mol for the PPDK reaction. If each compound is assigned only a fixed protonation level, based on the most abundant pseudoisomer at pH 7, the reaction looks as follows: PYR − + ATP 3− + P 2− i + H + PEP 2− + AMP − + PP 2− i which implies that the reaction energy increases with the pH, since the effect of one proton is RTln(10)·pH (Fig. 3). However, the PGC and PRC methods, which take the dissociation constants of the reactants into account, show that the response has a negative slope, which corresponds well to measured data (Fig. 3).

Comparing r G • estimations using the NIST-TECR database
In order to test the accuracy of the two methods described here, we use the NIST-TECR database as a benchmark. Given a measurement of K for a given reaction, we compare the predicted reaction energy (using any of the methods described in Section 3), to the observed reaction energy ( Fig. 1-equation I). The analysis of the PRC and PGC methods is given in Figure 4 and results in root mean squared error (RMSE) of 2.4 and 8.5 kJ/mol, respectively-equivalent to errors of about a factor of 3 and 30 in the estimation of K . Note that all of the different methods, those developed in the past as well as those presented here, have used some or all of the data that appear in NIST-TECR for training the group contributions or formation energies. As a result of these dependencies and the fact that the training procedures were not made public for the previously published methods, there is no way of performing an independent cross-validation for them. We found that when using the PRC method, there is a good fit between the regression data and the observed data in NIST-TECR (RMSE of 2.4 kJ/mol) while the methods based on group contributions do not fit the NIST-TECR data as well (RMSE of 8.5 kJ/mol for the PGC method and 9.9 kJ/mol for Jankowski et al. (2008)). The reason for this difference in accuracy could be attributed to the fact that group contribution is based on the simplifying assumption of independence between the contributions of groups to the overall f G • . In addition, there are more free variables in PRC than in PGC (407 compounds versus 99 groups). A summary of the analysis for the four different estimation methods based on the data provided by NIST-TECR is given in Table 1. The goodness of fit is given by the RMSE for each of the methods. These values should be compared with the baseline method which is to use the average r G • of each reaction in NIST-TECR across all its measurements. This baseline achieves the maximum accuracy but is limited in coverage to the scope of NIST-TECR. The baseline RMSE is 1.3 kJ/mol and is the average standard deviation of r G • per reaction.

Determining the reversibility of reactions
As an example of the usefulness of having a framework that makes all thermodynamic data available in one location and in an open format, we discuss the issue of reaction directionality that plays a pivotal and often problematic role in many metabolic models (Feist et al., 2007;Oberhardt et al., 2009) and has a crucial effect on their results. A reaction is called irreversible if its net flux flows in the same direction under all allowed physiological conditions. Some reactions are indisputably irreversible, for example the reaction of ribulose-bisphosphate oxygenation (promiscuously catalyzed by the enzyme RuBisCO) which has a r G • of -530 kJ/mol and therefore K = 10 93 . It is thus tempting to use a rule-of-thumb for determining whether a reaction is irreversible (Tanaka et al., 2003), by applying a threshold on its apparent equilibrium constante.g. K > 1000 (or K < 0.001 for irreversible reactions that always flow in the backward direction). This points to the fact that reactant concentrations are bound by physiological considerations and therefore a high-enough K would require too much of an imbalance in concentrations between substrates and products for reversing a reaction.
Aside from K itself, the reaction stoichiometry plays a crucial role in determining reversibility as well. For example, the fructose-bisphosphate aldolase reaction-fructose 1,6-P dihydroxyacetone-P + glyceraldehyde 3-P-has a r G • of 23 kJ/mol and a K of 9.23×10 −5 (Alberty, 2003). This might be considered irreversible (i.e. always flowing in the backward direction), but when the three reactants are at a concentration of 100 μM, a typical intracellular concentration for metabolites (Bennett et al., 2009), the r G is about 0.2 kJ/mol, very close to zero. The reason for the huge difference is due to the fact that this reaction has a different number of products and substrates, and in low-enough concentrations, the effect of the reaction quotient (Q ) on r G is significant. Therefore, any reversibility index should properly account for this stoichiometric effect. Furthermore, if a reaction involves many substrates and products, the dynamic range of r G • given the physiological metabolite concentrations can be much wider. That is, a reaction with one substrate and one product with r G • = 30 kJ/mol is much more irreversible than a reaction with three substrates, three products and the same r G • .
Previous studies (Feist et al., 2007) obtained binary reversibility annotations by checking if r G can attain positive and negative values at different physiological concentrations. Herein, we present a quantitative reversibility index,γ , which accounts for the effects of stoichiometry and physiological concentrations and defines a convenient metric for comparing the reversibility of reactions. The reversibility index is defined using the formulaγ ≡ K /Q 2/N where N is the total number of reactants (substrates plus products) and Q is the reaction quotient at characteristic physiological concentrations. For the purpose of calculating Q , all metabolites are taken to have a concentration of 100 μM (see Section 3 for details). The value ofγ signifies the fold change that all product and substrate concentrations must undergo in order to reverse a reaction. When considering a range ofγ <10 3 , which corresponds to allowing concentrations to span three orders of magnitude around 100 μM (∼3 μM-3mM), about 55% of the reactions are found to be reversible (see the Supplementary Material for a detailed statistical analysis).

DISCUSSION
The advances in metabolic modeling (Heinrich et al., 1991;Pfeiffer and Schuster, 2005;Oberhardt et al., 2009;Bar-Even et al., 2010), have created a need for accurate genome-wide values for reaction thermodynamic parameters. As metabolic network models for more organisms and cells emerge, it is increasingly important to have correct predictions for acidic and basic environments and for different ionic strengths. This requirement is most prominent when modeling organisms with multi-compartmented cells, and having to adjust the Gs to the conditions in each compartment. Most data that do exist are hard to access (e.g. in out-of-print books) and cover only a limited part of the scope of required compounds and reactions. Moreover, the attempts to use group contribution to expand this coverage are based on a closed implementation and are not extendable for outside users to add new group definitions or training examples.
In this article, we explored two approaches to estimating free energies. We show that the formation energies obtained using PRC achieve a good fit to the observed data (RMSE 2.4 kJ/mol), but do not provide genome-wide coverage of metabolic reactions (∼30% of E. coli model). As an alternative, PGC does provide free energy estimates for the majority (∼77%) of known biochemical reactions, but with larger errors (RMSE 8.5 kJ/mol). Therefore, a combination of the two methods might be beneficial, where PRC values are used whenever possible, and PGC is used to fill the gaps. There is, however, a potential problem with this approach as combinations of reactions can become inconsistent (e.g. stoichiometrically balanced cycles can have a non-zero G). The challenge of combining r G • estimations from different sources and estimation approaches in a unified and consistent manner requires an update to the methods described above (manuscript under preparation).
We hope that, with time, new measurements of reaction equilibrium constants will be published and used to improve the accuracy and coverage of both these methods (PGC and PRC). We thus join the plea of (Jankowski et al., 2008), who published a table of compounds that contain groups with yet unknown contributions.

CONCLUSION
We believe that the tools and data that enable thermodynamic analysis of biochemical systems should be easily and freely accessible. In addition to supplying r G • predictions as a table for use in metabolic models, we created a website (http://equilibrator.weizmann.ac.il) with a simple user interface that enables anyone to find reactions by chemical formula or enzyme name (Flamholz et al., 2011). The user can adapt the concentration of reactants and the conditions of the reaction.
The thermodynamics of biochemical reactions has a key role to play in our understanding and manipulation of metabolic pathways.
An integrated and open framework that combines accuracy and coverage will facilitate the wide use of this fundamental constraint by physics on the biochemistry of life.


extfieldsofstudy: Biology Medicine

text:
4-Methylenesterols from a Sponge Theonella swinhoei

Three new 4-methylenesterols, theonellasterol K (1), acetyltheonellasterol (2) and acetyldehydroconicasterol (3), along with two known sterols, theonellasterol (4) and theonellasterone (5), were isolated from the sponge Theonella swinhoei. The structures of these compounds were elucidated on the basis of their spectroscopic data and comparison of the NMR data with those of known analogues. Compound 1 exhibited significant cytotoxic activity against HCT-116, K562 and Molt 4 cancer cell lines.

Results and Discussion
The EtOAc extract of the freeze-dried specimen was fractionated by silica gel column chromatography and the eluted fractions were further separated utilizing normal phase HPLC to yield metabolites 1-5. The new compounds were given the trivial names theonellasterol K (1), acetyltheonellasterol (2) and acetyldehydroconicasterol (3). The known compounds were identified as theonellasterol (4) and theonellasterone (5). Theonellasterol K (1) was obtained as a white amorphous solid. The HRESIMS spectrum of 1 exhibited a pseudomolecular ion peak at m/z 481.3659 [M + Na] + , which established a molecular formula of C 30 H 50 O 3 , implying six degrees of unsaturation. IR absorptions were observed at 3362 and 3251 cm −1 , suggesting the presence of hydroxy groups in 1. The structure of this compound was deduced from its 13 C NMR and DEPT spectroscopic data (Table 1), which showed that the compound has 30 carbons, including six methyls, ten sp 3 methylenes, one sp 2 methylene, seven sp 3 methines (including one oxymethine), one sp 2 methine, three sp 3 quaternary carbons and two sp 2 quaternary carbons. From the 1 H NMR spectrum of 1, resonances of one olefinic methine proton (δ 5.68, d, J = 5.5 Hz), two olefinic methylene protons (δ 5.18 and 4.74, each s) and one oxygenated methine (δ 4.01, dd, J = 11.0, 5.0 Hz) were observed. Moreover, the 1 H NMR spectrum revealed the presence of one hydroperoxy proton resonating as a broad singlet at δ H 6.79. The planar structure and all of the 1 H and 13 C chemical shifts of 1 were elucidated by 2D NMR spectroscopic analysis, in particular 1 H-1 H COSY and HMBC experiments ( Figure 2). From the 1 H-1 H COSY correlations, it was found that one ring-juncture methine proton H-5 (δ 1.97) and one oxymethine proton H-3 (δ 4.01) exhibited allylic correlations with the exomethylene protons at C-30 (δ H 4.74 and 5.18). In addition, 1 H-1 H COSY spectral analysis established five partial structures of consecutive proton spin systems ( Figure 2). Further analysis of the HMBC correlations was employed successfully to establish the gross structure of 1 ( Figure 2). Thus, 1 was found to possess two double bonds at C-7/C-8 and C-4/C-30, one hydroxy group at C-3 and one hydroperoxy group at C-14. The relative configuration of 1, elucidated mainly from the NOESY spectrum, was compatible with that of 1 ascertained using molecular mechanics calculations (MM2), which suggested the most stable conformations, as shown in Figure 3. In the NOESY spectrum of 1, the NOE correlations between one of the methylene protons at C-11 (δ H 1.47) and both methyls (H 3 -18 and H 3 -19); H 3 -18 and H-20 as well as between H 3 -18, H 3 -19 and H-20 indicated that these protons adapt a β-orientation. H-5 was found to interact with H-3, but not with H 3 -19, revealing the β orientation of the hydroxy group at C-3. Moreover, the α-orientation of 14-OOH was further confirmed by the lower field chemical shift of H-17 (δ H 1.85, m). Furthermore, the chemical shifts of the side chain from C-20 to C-29 in 1 were nearly identical to those of 4 and 5. Thus, the structure of steroid 1 was established. After determining the structure of 1, we discovered that its molecular framework has been obtained as a known 4-methylenesterol conicasterol H, which was isolated previously from sponge Theonella swinhoei [18]. Acetyltheonellasterol (2) was isolated as a white powder with the molecular formula C 32 H 52 O 2 , which possesses seven degrees of unsaturation, as indicated by HRESIMS (m/z 491.3862, [M + Na] + ) and NMR spectroscopic data (Table 1). By comparison of the NMR data of 2 with those of 4, it was found that the 1 H and 13 C NMR data of 2 were very similar to those of 4, with the difference that 2 contains one more acetyl group relative to 4. The chemical shift of H-3 in 4 (δ H 4.01) was shifted downfield (δ H 5.15) in 2, suggesting that 2 is the 3-acetyl derivative of 4. We observed further that acetylation of 4 gave a product which was found to be identical to 2 by comparison of the physical and spectroscopic data. Thus, compound 2 was established as the 3-acetyl derivative of 4.
The new metabolite acetyldehydroconicasterol (3) was obtained as a white powder and possessed the molecular formula C 31 H 48 O 2 , as established from the HRESIMS and NMR data, implying eight degrees of unsaturation. Both the 1 H and 13 C NMR signals of 3 were found to be very closely related to those of compound 2, suggesting a very similar steroidal skeleton. By comparison of the NMR data of 3 with those of 2 (Table 1)

General Experimental Procedures
Optical rotation values were measured with a Jasco P-1010 digital polarimeter. IR spectra were recorded on a Varian Digilab FTS 1000 Fourier transform infrared spectrophotometer. The NMR spectra were recorded on a Varian Mercury Plus 400 FT-NMR (or Varian Unity INOVA 500 FT-NMR) instrument at 400 MHz (or 500 MHz) for 1 H-NMR and 100 MHz (or 125 MHz) for 13 C-NMR, respectively, in CDCl 3 . ESIMS were obtained with a Bruker APEX II mass spectrometer. Gravity column chromatography was performed on silica gel (230-400 mesh, Merck). TLC was carried out on precoated Kieselgel 60 F254 (0.2 mm, Merck) and spots were visualized by spraying with 10% H 2 SO 4 solution followed by heating. High-performance liquid chromatography (HPLC) was performed using a system comprised of a Hitachi L-7100 pump and a Rheodyne 7725 injection port. A preparative normal phase column (250 × 21.2 mm, 5 μm) was used for HPLC.

Animal Material
The specimen of Theonella swinhoei was collected by scuba divers at a depth of 15-20 m from coral reefs off the coast of Pingtung, Taiwan. A voucher specimen was deposited in the National Museum of Marine Biology and Aquarium, Taiwan (specimen no. 2011SP-2). Taxonomic identification was performed by Prof. Wen-Been Chang of the National Museum of Marine Biology & Aquarium, Pingtung, Taiwan.

Molecular Mechanics Calculations
Implementation of the MM2 force filed in Chem3D Pro software [24] was used to calculate the molecular models.


extfieldsofstudy: Chemistry Medicine

text:
Westphalen’s diol diacetate: 19(10→5)-abeo-5β-cholest-9-ene-3β,6β-diyl diacetate

The structure of the title steroid [alternative name: 3β,6β-diacetoxy-5β-methyl-19-norcholest-9(10)-ene], C31H50O4, confirms the generally accepted mechanism for the rearrangement of a cholestan-5α-ol derivative reported a century ago by Westphalen. The methyl group at position 10 of the starting material migrates to position 5 in the steroidal nucleus, while a Δ9 bond is formed, as indicated by the C=C bond length of 1.347 (4) Å. The methyl transposition leaves the 5R configuration unchanged, with the methyl oriented towards the β face. During the rearrangement, the steroidal B ring experiences a conformational distortion from chair to envelope with the C atom at position 6 as the flap. In the title structure, the isopropyl group of the side chain is disordered over two positions, with occupancies of 0.733 (10) and 0.267 (10). The carbonyl O atom in the acetyl group at C3 is also disordered with an occupancy ratio of 0.62 (4):0.38 (4).

The structure of the title steroid [alternative name: 3,6diacetoxy-5-methyl-19-norcholest-9(10)-ene], C 31 H 50 O 4 , confirms the generally accepted mechanism for the rearrangement of a cholestan-5-ol derivative reported a century ago by Westphalen. The methyl group at position 10 of the starting material migrates to position 5 in the steroidal nucleus, while a Á 9 bond is formed, as indicated by the C C bond length of 1.347 (4) Å . The methyl transposition leaves the 5R configuration unchanged, with the methyl oriented towards the face. During the rearrangement, the steroidal B ring experiences a conformational distortion from chair to envelope with the C atom at position 6 as the flap. In the title structure, the isopropyl group of the side chain is disordered over two positions, with occupancies of 0.733 (10) and 0.267 (10). The carbonyl O atom in the acetyl group at C3 is also disordered with an occupancy ratio of 0.62 (4):0.38 (4).

Related literature
For the initial report on the Westphalen rearrangement, see: Westphalen (1915). For applications in steroid synthesis, see: Rodig et al. (1961); Knights & Hanson (2004); Pinto et al. (2008Pinto et al. ( , 2009. For mechanistic aspects of this rearrangement, see: Kočovský & Č erný (1977); Kočovský et al. (1979); Kamernitskii et al. (1987). prepared under the Westphalen conditions in the cholestane series. Treatment with acetic anhydride and sulfuric acid promotes dehydration of the cholestan-5α-ol substrate, and the resulting carbocation is rearranged through the methyl 1,2-shift, before formation of the olefinic bond Δ 9 . The empty p orbital of the carbocation is electronically stabilized by the electron density of the O atom of the acetate at C6. With such a mechanism, a single epimer is expected as product, with the methyl at C5 oriented towards the β face. This stereochemistry has been confirmed in all the studied cases, but, surprisingly, the product of the very first report by Westphalen (1915) was never X-ray characterized.

Experimental
During a study on the optimization of reaction conditions for the Westphalen rearrangement carried out on 5-hydroxy-5α-cholestane-3β,6β-diyl diacetate, which was the substrate used by Westphalen, we obtained in 80% yield the rearrangement product (see Experimental). The molecular structure ( Fig. 1) and the absolute configuration for chiral centers are as expected. The saturated A and C rings have a chair conformation, while the B ring, which includes the C9═C10 double bond, is distorted to a half-chair conformation, close to an envelope with C6 as flap. Finally, the D ring adopts a half-chair conformation twisted on C13-C14. The observed conformation is indeed very similar to that described by Pinto et al. (2008) for the closely related derivative 3β-acetoxy-6β-hydroxy-5β-methyl-19-norcholest-9(10)ene. An overlay between the steroidal nucleus of this structure and that of the title molecule gives a r.m.s. deviation limited to 0.156 Å, the largest deviation arising from the side chain. This group presents a degree of flexibility, as reflected by the disorder detected in the title compound for the isopropyl group. The same conformation for the A-D ring system was also observed for a Westphalen product in the pregnane series (Pinto et al., 2009).
Experimental 5-Hydroxy-5α-cholestane-3β,6β-diyl diacetate (6.1 g, 12.09 mmol) was treated with acetic anhydride (150 ml) for 20 min at 363 K. Then, NaHSO 4 was added (1.74 g, 12.78 mmol) and stirred until complete consumption of the starting material (ca. 30 min). The mixture was poured in a 500 ml Erlenmeyer with pyridine and ice and vigorously stirred. The product, which deposited on the glass vessel, was dissolved in ethyl acetate, and then washed with saline solution, 5% HCl, distilled water, and finally 10% NaHCO 3 . This phase was dried over Na 2 SO 4 , and the solvent evaporated under reduced pressure. The crude (80% yield) was chromatographed over silicagel with petroleum ether and ethyl acetate (95:5) and crystallized from ethyl acetate. When KHSO 4 was used, the Westphalen compound was obtained in 35% yield.

Refinement
The isopropyl group in the lateral chain was found to be disordered over two positions, C251/C261/C271 and C252/C262/C272, and occupancies converged to 0.733 (10) and 0.267 (10), respectively. Both parts were restrained to have similar displacement parameters and geometry (SIMU and SAME restraints; Sheldrick, 2008). On the other side of the molecule, the carbonyl O atom in the acetyl group at C3 is also disordered over two sites, O301 and O302, with occupancies 0.62 (4) and 0.38 (4). All H atoms were placed in calculated positions and refined as riding to their carrier C atoms. C-H bond lengths were set to 0.96 (methyl) 0.97 (methylene) or 0.98 Å (methine), and isotropic parameters for H atoms were calculated as U iso (H) = xU eq (carrier C) with x = 1.5 (methyl) or x = 1.2 (methylene and methine groups).
Friedel pairs (714) were merged, and absolute configuration assigned by fixing the configuration of known chiral centers.

Figure 1
ORTEP view of the title molecule, with displacement ellipsoids for non-H atoms at the 30% probability level. For disordered parts, sites with occupancy less than 0.5 have been omitted.