Dynamic approach to predict pH profiles of biologically relevant buffers.

Recently, dynamic approach has been applied to determine the steady state concentrations of multiple ionic species present in complex buffers at equilibrium. Here, we have used the dynamic approach to explicitly model the pH profiles of biologically relevant phosphate buffer and universal buffer (a mixture of three tri-protic acids such as citric acid, boric acid and phosphoric acid). The results from dynamic approach are identical to that of the conventional algebraic approach, but with an added advantage that the dynamic approach, allow for the modelling of complex buffer systems relatively easy compared to that of algebraic method.

Introduction

Buffer preparation is an integral part of several in-vitro and in-vivo experiments performed in biological studies. Buffer by definition is a mixture of weak acid and its conjugate base, or a mixture of weak base and its conjugate acid. The pK of the weak acid or the pK of the weak base defines the buffering capacity of the buffers used for in the experiments. Simple buffers are prepared by titrating specified quantity of acid/base (mono-, di- or tri- acidic/alkaline) against alkali/acid till the desired pH is achieved. It is necessary that the final pH should be within the buffering capacity of the acid/base being considered, for the buffer to be effective. Biologically relevant buffers such as phosphate buffers are prepared by titrating equi-molar concentration of mono hydrogen phosphate against dihydrogen phosphate or vice-versa till a desired pH is reached. Complex buffer such as universal buffer are prepared by titrating a mixture of two or three tri-protic acids with an alkali or tri-sodium salt of phosphoric acid [1].

Two mathematical approaches namely, the algebraic [2], [3] and the dynamic approach [4], [5], [6] are available to predict the pH profiles of complex buffer systems. In this work, we explicitly derive the analytical expression for phosphate and universal buffer using algebraic as well as differential dynamic approach. The derived models were used to back-predict the pH profiles of phosphate and universal buffer solutions reported earlier [7]. The theoretical predictions are in good agreement with the experimental data, validating the accuracy of the proposed models.

Theory

Consider the volumetric titration of a weak mono-protic acid () in the presence of strong alkali such as . At the outset, the weak acid exists in the form of both unionized form and its ionized form . The equilibrium between these two states can be written as (Eq. (1)),

, are the kinetic rates for the dissociation and association of the weak acid, respectively. The irreversible dissociation of alkali during the titration is given by (Eq. (2)),

is the kinetic rate for the irreversible dissociation of the alkali. Unlike, weak acid, the dissociation of the strong base/acid are represented through irreversible reactions. Finally the dissociation of water is given by (Eq. (3)),

, are the forward and reverse kinetic rates for the dissociation of water, respectively. The pH profile of the weak acid upon addition of alkali results in a typical weak acid-base titration curve containing buffering region (where, pH=pKa) and equivalence point [8], [2], [3]. To predict the pH profile of a titration, an analytical equation that relates pH with the concentration of the alkali is necessary [9]. The derivation for such titrations are available in several literatures, text-books and online resources [10] (please refer supporting information). The first step in the derivation involves, expression of the mole fractions of [HA] and [A] in terms of total acid [A] and the equilibrium constant K, using law of mass balance equation [HA]+[A]=[A]. The second step involves, assumption of charge neutrality of the solution, i.e, [H]+[Na] =[OH]+[A]. Re-expression of all variables in the charge balance equation, in terms of [H] would result in a polynomial equation. Based on the order of the polynomial equation, multiple roots are available. Choosing the correct root value that is physically relevant would yield the [H] ion concentration. The negative logarithm to the base 10 of [H] is defined as the pH value (pH=-log [H]). By varying [Na], the titrated alkali concentration, different polynomials are obtained. Each polynomial equation is solved separately to determining the pH profile of the titration curve.

Alternatively, pH can also be predicted using dynamic approach, where we first assume a kinetic model and then frame a set of differential equation based on the model [6] (please refer Supporting information).

In the above differential equations from (4), (5), (6), (7); [HA] and [A] are the unionized and the ionized species of the weak acid, respectively. [H] is the proton generated by dissociation of acid as well as water. [OH] is the hydroxyl ion generated from dissociation of alkali as well as water. [NaOH] is the alkali component (sodium hydroxide) used in the titration. Since, [NaOH] dissociates completely to [Na] and [OH], immediately after it's addition to the solution, the dissociated [OH] arising from [NaOH] can be directly added to the initial concentration of [OH], present at the onset of the titration. Thereby, the rate equation for [NaOH] can be explicitly excluded from the above model.

In order to obtain a steady-state concentration, the above set of coupled differential equations are integration over a long period of time, such that, the equilibrium state had been reached. The [H] concentration obtained after equilibrium is achieved can be used to calculate the pH profile of the titration. The simulations of mono-protic, di-protic, tri-protic, phosphate buffer and universal buffer titrated with NaOH is shown in Fig. 1. Additionally, the titration of amino-hydroxyl compound with HCl is also presented in Fig. 1. The simulations carried out using algebraic as well as dynamic approach generated identical results. For the simulation of mono-protic acid, the rate constants , ,, , were set to 10−4, 1.0, 10–14, 1.0 s−1, respectively. The matlab [11] codes for simulation of mono-protic acid is provided in the supporting information. The simulation of di-protic acid was carried out using rate constants, , , , ,, , whose values set to 10−4, 1.0, 10−9, 1.0, 10–14,1.0 s−1, respectively. The , are the , rate constants for the second ionisable moiety of the di-protic acid. The simulation of tri-protic acid was carried out using rate constants, , , , , , ,, , whose values set to 10−3, 1.0, 10−6, 1.0, 10−9, 1.0, 10–14,1.0 s−1, respectively. The , are the , rate constants for the third ionisable moiety of the tri-protic acid. The simulation of mono-hydroxyl alkali was carried out using rate constants, , ,, , whose values set to 10−9.24, 1.0, 10–14,1.0 s−1, respectively. The , are the , rate constants for the dissociation of alkali component. The matlab codes for simulation of mono-alkaline species is provided in the supporting information.

Fig. 1

The simulation of pH profile obtained for titration of mono-protic (A), di-protic (B), tri-protic (C), mono-alkaline (D), Phosphate buffer (E), Universal buffer (F). NaOH was considered as a titrant for (A), (B), (C), (F). HCl was considered as titrant for (D). Na2HPO4 was considered as titrant for (E). For simulation, pKa was set to 4 for mono-protic acid, 4,9 for di-protic acid, 3,6,9 for tri-protic acid, 4.6 or pKb=9.245 for mono-alkaline, 2.12,7.21,12.32 for phosphate buffer. For universal buffer the pKa values were set to 2.12, 7.21, 12.32 for phosphoric acid, 3.08, 4.72, 5.4 for citric acid, 9.23, 12.4, 13.3 for Boric acid. The simulation of phosphate buffer assumes titration of NaH2PO4 with incremental volume of Na2HPO4. Both the algebraic and dynamic approach yielded identical profiles for all the simulations.

Phosphate buffer

Kinetic model for phosphate buffer

Phosphoric acid is a tri-protic acid with three pK values. Preparation of phosphate buffer involves titration of equi-molar concentration of mono hydrogen phosphate into di-hydrogen phosphate or vice versa. We assume the following kinetic model to predict the pH profile of the phosphate buffer,

In the above Eqs., (8), (9), (10), [HPO], [HPO], [PO], represent the mono-, di-, trivalent phosphate species generated by dissociation of the phosphoric acid [HPO].

In the above Eqs. (11), (12), (13), is the acid component added in the form of mono-sodium salt of phosphoric acid. is the alkali component added in the form of di-sodium salt of phosphoric acid. [H] is the proton generated during the dissociation of acid as well as water; and [OH] is the hydroxyl ion generated by dissociation of alkali as well as water. , , , are the forward kinetic rates and , , , are the reverse kinetic rates for the equilibriums in Eqs. (8), (9), (10), (13), respectively. , , in Eqs. (11), (12) are the irreversible dissociation rates of , , respectively. The simulation of phosphate buffer was carried out using rate constants, , ,, , , ,, , whose values set to 10−2.12, 1.0, 10−7.21, 1.0, 10–12.32, 1.0, 10–14,1.0 s−1, respectively. The irreversible dissociation rates of both mono and di-sodium phosphate salts, , were set to 102 s−1 (Fig. 1).

Algebraic method for phosphate buffer

The charge neutrality equation for such a system can be written as (Eq. (14)),

Since, and , the mono-, di-sodium salt of phosphate, dissociate completely, its concentration can be directly substituted in to the above equation. Further, and , represent the acidic and alkali components of the titration. By re-expressing all the other components such as [ OH], [HPO], [HPO], [PO], in terms of [H], the above equation simplifies to a polynomial of order 5 in [H] (Eq. (15)).

, , and are the equilibrium constants obtained from the equilibrium conditions, represented by Eqs., (8), (9), (10), (13), respectively. The ratio of forward kinetic rate to reverse kinetic rate yields the equilibrium constants; i.e., ,, . By solving the polynomial, Eq. (15), for a given and , the concentration of and , and selecting a physically relevant root, the pH can be determined. By varying A and B, in Eq. (15), the pH profile of a complete titration can be obtained. (Matlab code to symbolically derive Eq. (15) is available in the supporting information).

Dynamic approach for phosphate buffer

To predict the pH for a titration curve using differential approach, we frame the following system of differential equations based on the kinetic model (Eqs. (8), (9), (10), (11), (12), (13)),

The contribution of the dissociation of mono-sodium salt, [NaHPO]0 and di-sodium salt, [NaHPO] are added directly to the initial concentrations of [HPO] and [HPO], respectively. The above set of coupled differential equations (Eqs.(16), (17), (18), (19), (20), (21)) can be numerically integrated for a given initial condition/concentration set at 0, [NaHPO], [NaHPO], 0, 10−7, 10−7, respectively. [NaHPO], [NaHPO], are the initial concentrations of mono-sodium phosphate and di-sodium phosphate salts, respectively, which constitutes the acidic and alkali components present in the mixture. The concentrations of [NaHPO] and [NaHPO], were varied for each simulation to obtain the entire pH profile of the titration data. The simulation of phosphate buffer using algebraic and dynamic approach is shown in Fig. 1(E).

Universal buffer

Kinetic model for universal buffer

Universal buffer contains a mixture of three tri-protic acids namely, phosphoric acid, citric acid and boric acid. The preparation of universal buffer involves, preparing a mixture of citric and boric acid at appropriate concentrations and titrating it against tri-sodium salt of phosphoric acid to a desired pH value.

We assumed the following kinetic model for the universal buffer:

In the above Eqs. (22), (23), (24), (25), (26), (27), (28), (29), (30), (31), (32), [HPO], [HPO], [PO], are the mono-, di-, tri-, anionic species generated from phosphoric acid [HPO], respectively. Na3PO4 is the alkali component added in the form of tri-sodium salt of phosphoric acid. [HC], [HC], [C], are the mono-, di- and tri-, anionic species generated from citric acid [HC], respectively. [HB], [HB], [B], are the mono-, di- and tri-, anionic species generated from boric acid [HB], respectively. [H], is the proton generated by dissociation of acids as well as water and [OH], is the hydroxyl ion generated by dissociation of alkali as well as water. Except for Eq. (25), which is irreversible and leads to complete dissociation of [NaPO] with rate constant, , in solution, all the other reactions are reversible by nature. , , , , , , , , ,, are the forward kinetic rates and , , , , , , , , ,, are the reverse kinetic rates for the equilibriums in Eqs. (25), (27), (28), (30)–(33), (35), (36), respectively. The simulation of universal buffer was carried out using rate constants, , , , , , , , , , , , , , , , , ,, , , , whose values set to 10−2.12, 1.0, 10−7.21, 1.0, 10–12.32, 1.0, 102 ,10−3.08, 1.0, 10−4.74, 1.0, 10−5.4, 1.0, 10−9.23, 1.0, 10–12.4, 1.0, 10–13.3, 1.0, 10–14, 1.0 s−1, respectively. The rate constants through , corresponds to kinetic dissociation of phosphoric acid, through , corresponds to kinetic dissociation of citric acid, through , corresponds to kinetic dissociation of boric acid. The rate constant, , represents the irreversible dissociation of tri-sodium phosphate, which was considered to be the sole alkali component of the titration (Fig. 1).

Alternatively, universal buffer can also be prepared by titrating a mixture of phosphoric, citric and boric acid with NaOH. The kinetic model for such a system will differ by inclusion of dissociation of [NaOH] instead of [NaPO]. The simulation of universal buffer system with NaOH as alkali titrant is shown in Fig. 1(F). In the present work, we will be considering universal buffer which is prepared by titrating mixture of citric and boric acid with NaPO, as the sole alkali component.

Algebraic model used for universal buffer

The charge neutrality equation for the above kinetic model is given by Eq. (33),

By re-expressing all the components such as [OH], [HPO], [HPO], [PO], [HC], [HC], [C], [HB], [HB], [B] in terms of [H+] the above equation simplifies to a polynomial equation of order 11 expressed in [H]. The explicit derivation using matlab symbolic functions is available in the supporting information. By solving the polynomial equation and selecting the physically relevant root, the pH can be determined.

Dynamic approach for universal buffer

Dynamic approach requires the following set of differential equations based on the above kinetic model,

Since, [NaPO] the only alkali component present in the system, dissociates completely into solution as phosphate ion, the contribution of the dissociation of [NaPO] can be directly added to the initial concentration of [PO]. During numerical integration, the initial concentrations/conditions for each species starting from Eqs. (34), (35), (36), (37), (38), (39), (40), (41), (42), (43), (44), (45), (46), (47) are set to 0, 0, 0, [NaPO], [HC], 0, 0, 0, [HB], 0, 0, 0, 10−7, 10−7 respectively. [NaPO], [HC], [HB], are the total concentrations of tri-sodium phosphate, citric acid and boric acid, respectively, present in the mixture. The concentration of [NaPO], was varied for each experiment to simulate the experimental data.

Results

Experimental pH profiles for the standard phosphate buffer within the pH range 6–7.2 (Table 1) and 5.8–8 (Table 2) are available in molecular biology protocols [7]. In this experimental protocol, 0.1/1.0 M solution of mono sodium phosphate was mixed with 0.1/1.0 M di-sodium phosphate at specified ratios to prepare buffer solutions of desired pH [7]. At the outset, we tried to simulate the pH profile that would closely resemble the experimental pH profile, using the standard pK of phosphoric acid, i.e. 2.1, 7.2, 12.0. The simulated profile and the experimental profile showed discrepancy due to the fact that the standard values obtained from literature need not necessarily correlate with the optimized set of pK values that pertains to the experimental condition. In order to minimize the difference between the experimental and simulation we directly fitted the experimental data to the phosphate model by varying the pK values. In-house written matlab codes were used to perform the optimization of pK values. The curve fitting yielded an optimized pK values of 2.239, 6.864, 11.627 for phosphoric acid in a phosphate buffer prepared using 0.1 M stock solution. Similarly the optimized, pK values of 2.15, 6.8537 and 8.768 were obtained for phosphoric acid in a phosphate buffer prepared using 1 M stock solution. The pH profiles optimized using dynamic approach were in excellent agreement with the experimental values with a maximum deviation of 0.001 pH units at pH 7.6 and 0.015 pH units at pH 7.2 for 0.1 M and 1 M buffer solution, respectively (Fig. 2).

Table 1

The pH profile of the 1 M phosphate buffer prepared by mixing different ratios of 1 M NaH2PO4 and 1 M Na2HPO4, The experimental values are model fitted using dynamic approach and the back-calculated pH values based on optimized pKa values are tabulated below.

Solution "A" (mL)	Solution "B" (mL)	1 M Phosphate buffer (pH)
(1 M NaH2PO4)	(1 M Na2HPO4)	Experimental	Back calculated
877	123	6.0	5.989
850	150	6.1	6.088
815	185	6.2	6.196
775	225	6.3	6.302
735	265	6.4	6.395
685	315	6.5	6.498
625	375	6.6	6.611
565	435	6.7	6.716
510	490	6.8	6.807
450	550	6.9	6.905
390	610	7.0	7.003
330	670	7.1	7.102
280	720	7.2	7.188

Table 2

The pH profile of the 0.1 M phosphate buffer prepared by mixing different ratios of 0.1 M NaH2PO4 and 0.1 M Na2HPO4, The experimental values are model fitted using dynamic approach and the back-calculated pH values based on optimized pKa values are tabulated below.

Solution "A" (mL)	Solution "B" (mL)	0.1 M Phosphate buffer (pH)
( 0.1 M NaH2PO4)	( 0.1 M Na2HPO4)	Experimental	Back calculated
6.80	93.20	8.00	8.000
10.40	89.60	7.80	7.799
15.50	84.50	7.60	7.601
22.60	77.40	7.40	7.399
31.60	68.40	7.20	7.200
42.30	57.70	7.00	7.000
53.70	46.30	6.80	6.800
64.80	35.20	6.60	6.600
74.50	25.50	6.40	6.399
82.20	17.80	6.20	6.201
88.00	12.00	6.00	6.000
92.10	7.90	5.80	5.800

Fig. 2

The experimental (•) and model fitted (—) pH profile of phosphate buffer prepared by mixing different volumes of equimolar concentrations of NaH2PO4 and Na2HPO4. (A) Shows the pH profile of 0.1 M phosphate buffer prepared by mixing mono and di sodium salt of phosphate. The experimental data was fitted to phosphate buffer model and optimized pKa values such as 2.239, 6.864, and 11.627 were obtained. (B) Shows the pH profile of 1 M phosphate buffer prepared by mixing mono and di sodium salt of phosphate. The experimental data was fitted to phosphate buffer model and optimized pKa values such as 2.15, 6.8537 and 8.768 were obtained.

Titrated pH profile for universal buffer within the pH range 2–12 (Table 3) is obtained from earlier literature [1]. As seen for phosphate buffer, the simulated pH profile of universal buffer also showed discrepancy with the experimental value when standard pK values (Citric acid: 3.1, 4.7, 5.4; Boric acid: 9.2, 12.4, 13.3) were used, owing to the fact that pK values are sensitive to ionic concentration and temperature [1]. To minimize the difference between the simulated profile and the experimental data, we optimized the pK values through curve fitting using dynamic model. The optimized pK values for phosphoric acid were 3.259, 7.029, 12.233. The optimized pK values for citric acid were, 2.628, 4.431, 5.620, and the optimized pK values for boric acid were and 9.770, 12.507, 12.92. The fitted pH profile using optimized pK values for phosphoric, citric and boric acid agreed well with the experimental data with a maximum deviation of only 0.22 pH unit for pH 10.5 (Fig. 3).

Table 3

The pH profile of the universal buffer prepared by mixing different ratios of 0.2 M Boric Acid+0.05 M Citric Acid and 0.1 M Na₃PO4, The experimental values are model fitted using dynamic approach and the back-calculated pH values based on optimized pKa values are tabulated below.

Solution "A" (mL)	Solution "B" (mL)	Universal buffer (pH)
(0.2 M Boric Acid+0.05 M Citric Acid)	0.1 M Na3PO4	Experimental	Back calculated
195	5	2.0	1.918
184	16	2.5	2.607
176	24	3.0	2.988
166	34	3.5	3.475
155	45	4.0	4.010
144	56	4.5	4.510
134	66	5.0	5.000
126	74	5.5	5.449
118	82	6.0	5.963
109	91	6.5	6.580
99	101	7.0	7.101
92	108	7.5	7.441
85	115	8.0	7.874
78	122	8.5	8.540
69	131	9.0	9.181
60	140	9.5	9.622
54	146	10.0	9.933
49	151	10.5	10.279
44	156	11.0	10.809
33	167	11.5	11.597
17	183	12.0	12.047

Fig. 3

The experimental (•) and model fitted (—) pH profile of universal buffer prepared by titrating Na3PO4 (0.1 M) against a mixture of citric acid (0.05 M) and boric acid (0.2 M). The optimized pKa values for phosphoric acid were 3.259, 7.029, 12.233, citric acid were 2.628, 4.431, 5.620 and boric acid were 9.770, 12.507, 12.920.

Discussion

It is natural for the simulated pH profiles to deviate from the experimental pH profiles, when standard pK values are used in the calculation. Appropriate correction factors are to be incorporated to the standard pK values before any comparison is made between the simulation and the experimental pH profiles. Such discrepancies arise largely due to the dependence of pK on physical factors such as ionic strength and temperature. The explicit expression for pK is given as pK=–log(K); where K is the equilibrium constant of the dissociation process [12]. In general, the K is calculated from the ratio of concentrations of product to reactants, i.e. . Here, P are the products of the reaction, R are the reactants of the reaction; n and m are the number of reactants and products involved in the reaction, respectively. , represents the multiplication of each terms whose index goes from i=1 to m or j=1 to n for reactant and products, respectively. Practically, due to ionic nature of the P and R, the concentrations [P] and [R], ought to be replaced with activities and [13], which are represented as the multiplication of activity coefficient (γ) with the respective concentrations ([P]) i.e. = γ [P] and = γ [R]. To obtain accurate pK values which compares well with the experimental results for the simulations, , the thermodynamic equilibrium constant, should be used instead of K. The mathematical expression that relates with K is, . If , the cumulative activity coefficient, can be determined experimentally, then the correction factor due to ionic strength at high concentrations can be easily incorporated in to the pK value. In the case of phosphate buffer, the standard pKa value for the dissociation of di-anionic phosphate to tri-anionic phosphate is around 12.0; whereas, the optimized pKa value for 0.1 M and 1 M phosphate buffer were, 11.6 and 8.8, respectively. A significant decrease of ~2.8 pH units observed for the optimized pKa value in the case of 1 M phosphate buffer, in comparison with the standard value, is largely due to the increased ionic effect seen at high concentration of the buffer solutions.

In addition to ionic activity, temperature also affects the value of thermodynamic equilibrium constant, [14]. The relationship between temperature T and is given by [10]. Where, , is the Gibbs free energy, R is the Gas constant. Temperature dependent correction factors are also needed be incorporated in to pK, when and T are known. By making appropriate corrections to pK based on ionic activity and temperature effect, pH profiles that closely resembles the experimental values can be realized. In our case, we directly fitted the experimental data to the buffer models based on dynamic approach to obtain optimized pK values, circumventing the need for calculating the correction factors due to ionic and temperature effect. The optimized pK values are infact a self-corrective pK that implicitly takes into account the combined effects of ionic strength and temperature.

During numerical simulation of pH profile using dynamic approach, in addition to using optimized pK values, it is also necessary to evolve the system of differential equations for a longer period of time. If we let the differential equations to evolve for a shorter time period, the resulting solutions will only represent a pre-steady state condition which is still dynamically varying with time rather than the final equilibrium state. Hence, it is quint essential to evolve the system of differential equations for a longer period of time (e.g.>1015 s) so as to attain accurate pH values at equilibrium condition.

Conflict of interest

The authors declare no conflict of interest.