ARS Home » Southeast Area » Raleigh, North Carolina » Food Science and Market Quality and Handling Research Unit » Research » Software » Download

Buffer Capacity Model

Buffer capacity modeling software: description of methods and software available for download.

Fred Breidt

fred.breidt@usda.gov

USDA/ARS Microbiologist

Introduction

Buffer capacity modeling methods and software have been developed (Price et al., 2020). This technology was applied to estimate the effects of low acid salad dressing ingredients on dressing pH (Longtin et al., 2020). The models may be useful for estimating the effects of ingredients on pH changes in a variety of acid or acidified foods, as well as for food safety and product development purposes. Software is available for download at link below. The software was developed using Matlab programming environment (version 2019b) and requires the Matlab Optimzation toolbox. Please address all questions to Dr. Fred Breidt at: fred.breidt@usda.gov.

Predicting pH from a matrix of buffer concentrations and pKs

To determine pH from a mixture of monoprotic weak acids and bases, the molar concentrations Ca and Cb for the acids and bases and the corresponding equilibrium constants Ka or Kb values are required. An algebraic equation may be derived from the combined weak acid equilibrium and charge balance equations to solve for [H+] (and therefore pH) with Kw being the equilibrium constant of water (10-14), as described by Butler & Cogley (1998, pp. 123):

acid term                   base term           hydroxyl [OH-] term
1) [H+] = (CaKa/(Ka + [H+])) - (Cb[H+]/([H+] + Kb)) + Kw/[H+]

To solve Equation 1 for [H+], the solution of a 4th order polynomial is required. Because the order of the polynomial equation increases with additional acid or base terms, numerical solutions may be used for estimating pH. Generalizing for multiple monoprotic acids and bases as well as rearranging and adding terms for the concentration of the salts of acids or bases Ani, (a positive value for anions) or Cti (negative value for cations) results in:

2)   0 = ∑(CaiKak/(Kai + [H+])) - ∑(Cbi[H+]/([H+] + Kbi)) + Kw/[H+] - [H+] + ∑ (Ani + Cti)

An algorithm was developed using Matlab (CalcpH_AB.m, Table 1) to obtain a numerical solution to Equation 2 for a given set of concentration, pK and An or Ct values using Newton’s minimization method as suggested by Butler & Cogley (1998). A defined tolerance (10-10 for the function value, or 1000 iterations) was used as the stopping criterion for the minimization algorithm. For pH estimation, the equilibrium constants for acids and bases (transformed by the negative log10 to pK values) were adjusted as needed for ionic strength due to added sodium calcium salts using the Davies equation (AdjpK_AB.m, Table 1) as suggested by Butler & Cogley (1998, pp. 94).

Titrations to generate BC curves

To generate buffer capacity curves, we used an automatic titrator with dual cylinders, one for NaOH, and one for HCl (Model 902, Hanna Instruments, Smithfield, RI, USA). The dynamic dosing feature was used to generate smooth titration curves from the initial solution pH up to pH 12 with NaOH, or down to pH 2 with HCl. Acid and base concentrations depended on the solution being titrated. Data for the acid and base titration curves of 50 ml solutions (initial volume) were exported from the titrator as a text file and processed using a custom Python script (GetCurve.py, F. Breidt, unpublished) to generate a comma delimited data files containing two columns: the volume of acid or base added and the resulting pH. The data were then imported into Matlab with the built in csvread.m function and saved as a Matlab workspace variable. Three additional workspace variables are needed: the concentration of the acid and base used for titration and the initial volume for the titration. A parameter file (text file, example shown in Table 2),  Matlab version R2019a was used for processing the data. A single Matlab “Live-Script” file was used to sequentially call the functions and process the data (example file: DriveFit_LiveSheet.mlx, Table 1). Model parameters are listed in Table 2. Parameters were stored in a text file and were imported as workspace variables at the start of the LiveSheet file (paramfile.txt, Table 1). To generate a BC curve from the titration data, an iterative step-wise derivative was calculated (Butler & Cogley, 1998, pp. 133) with the acid or base concentration used and volume titrated (function Data2Beta.m, Table 1):

3)    β = ∆(Acid or Base)/∆pH

where ∆Acid or ∆Base represented the change in normality of acid or base in the solution being titrated (based on the volume added at each titration step and the concentration of the HCl or NaOH) with the resulting pH change ∆pH. A minimum change (typically 0.02 pH units) in pH was defined to prevent deviations in the derivative if an air bubble was in the acid or base liquid delivery line of the titrator, or some other error occurred during titration (minDeltapH in Table 2). Prior to further processing, BC curves were manually trimmed by setting trim_beg or trim_end (Table 2) as needed, so the ends of each BC curve, typically at pH 2 and pH 12, had similar BC values. In addition, gaps in the BC greater than a set pH value, typically 0.3 pH units or greater (minGap, Table 2), were automatically filled with BC data points at a specified interval of 0.1 pH units (increment, Table 2) using a linear function between the end points of the gap (fillgap.m, Table 1). These steps were necessary to facilitate the subsequent use of curve fitting algorithms.

Modeling BC curves

Curve fitting for BC data was accomplished using a trigonometric least squares regression method (Eubank, 1990; Newbery, 1970), based on the equation:

4)    F(x) = Bo + A1sin(αx) + B1cos(αx) + A2sin(2αx) + B2cos(2αx) + A3sin(3αx) ...

where F(x) represented the BC value for a given pH (represented as the variable x), Ai and Bi were scalar parameters, and α was a multiplier (0.5) used for x to fit pH values between pH 2 and pH 12. The model was fit to the BC data by minimizing the error sum of squares:

5)    SS(E) = ∑(F(x) - Y(x))2

where F(x) was the trigonometric series model (Equation 4), and Y(x) was the BC curve generated from Equation 3. To estimate parameters for the trigonometric model, a square matrix (M) of partial derivatives of F(x) with respect to each parameter (Pi) of Equation 5 was generated:

6)   ∂SS(E)/∂(Pi) = 0

An optimized parameter vector (Po = B0, A1, B1, A2, B2, …) for the terms in Equation 4 was then obtained using Matlab’s ‘mldivide’ operator to solve the matrix equation M*Po = Yi, where the vector Yi consisted of the corresponding Y(x) term for each partial derivative equation (implemented as SCBCfit.m, Table 1). Typically, 15 Ai and Bi parameters (defined as the parameter “Order”, Table 2) were used for the solution. The resulting trigonometric series model with optimized parameters (Po) represented a mathematically defined smoothed curve representing the BC data and was used as a template for modeling the BC curve with the BC model as described by Butler & Cogley (1998, pp134):

7)    β = 2.303*(∑(CiKi[H+]/([H+] + Ki)2) + Kw/[H+] + [H+])

where Ci and Ki are the concentration and dissociation constants for weak acids or bases that contribute to buffering in the solution that was originally titrated and Kw is the equilibrium constant for water.

A matrix of Ci and Ki parameter values were obtained for the optimized fit of Equation 7 to the trigonometric function, using the constrained nonlinear minimization algorithm fmincon.m (Optimization Toolbox, Matlab) to minimize squared error values (implemented as SimplexBCPK_DF.m, Table 1). Starting values for the minimization algorithm were chosen by using a predefined number of pK values, typically 7 (parameter NpKs, Table 2), evenly distributed between the minimum and maximum pH values for the titration curve with the corresponding concentrations estimated (assuming pH = pK) from BC value from the trigonometric regression model. The NpKs parameter value could also be adjusted manually based on the complexity of the BC curve. Concentrations were constrained to positive values only, and Ki values (transformed to pK values) were constrained between the upper and lower pH values for the BC curve, typically pH 12 to pH 2 (defined as UB and LB, respectively, in the parameter file Table 2). For the resulting concentration-pK matrix, all pK values that were within a set tolerance of each other, typically 0.2 pH units (pK_tol, Table 2), were combined and the corresponding concentration values (Ci values) were added. The result was an Nx2 matrix of concentration and pK values for each ingredient. The estimated pH from the matrix was then determined using Equation 2 with An and Ct equal to zero, with corrections of the pK due to the NaCl concentration (NaClpercent, Table 2). An additional parameter (WaterSalt, Table 2) was used for correcting the BC of water for the final BC graph. By default, pK values for buffers identified by the model that were greater than pH 7 were considered bases for pH estimation (by Equation 2), and pK values for buffers with a pK less than or equal to pH 7 were considered acids.

To model salts of acids or bases for pH prediction from the BC model data, a positive or negative value representing the sum of An and Ct (a positive anion or negative cation value) was used to minimize the squared difference between the observed and predicted pH using Equation 2 (GetAdjC.m, Table 1). The magnitude of this value was essentially the millimolar concentration of the salt of the acid or base. If no acid or base salts were present, the magnitude of this “ion” value represented the error in the model for the predicted BC curve. If known acid or base salts were present, this value could be set in the parameter file (parameter adjC, Table 2).

References:

Butler, J. N., & Cogley, D. R. (1998). Ionic equilibrium: Solubility and pH calculations (pp. 94-201). New York, NY: John Wiley and Sons.

Eubank, R.L. & Speckman P. (1990). Curve fitting by polynomial-trigonomietric regression. Biometrika 77(1):1-9.

Newbery, A. C. R. (1970). Trigonometric interpolation and curve fitting. Mathematics of Computation, 24(112), 869-876.

Longtin M, Price RE, Mishra R, Breidt F. 2020. Modeling the buffer capacity of ingredients in salad dressing products. J Food Sci 85(4):910-917.

Price RE, Longtin M, Conley Payton S, Osborne JA, Johanningsmeier SD, Bitzer D, Breidt F. 2020. Modeling buffer capacity and pH in acid and acidified foods. J Food Sci 85(4):918-925.

Table 1. Matlab program dependencies for the model.

 Program name Description DriveFit_LiveSheet.mlx Matlab "LiveSheet" script file to sequentiall call processing algorithms Data2Beta.m Process titrator data to generate BC curve fillgaps.m Fill gaps in BC curves with a given increment and linear model SCBCfit.m Implement the trigonometric regression algorithm SimplexBCPK_DF.m Fit the buffer model to the trigonometric regression curve Beta2Conc.m Convert a BC value to a molar concentration for a buffer BetaModel_AB.m Calculate a BC curve from a matrix of buffer concentration and pK value CalcpH_AB.m Calculate pH from a matrix of buffer concentration and pK values ConbineABs.m Merge buffers with closely related pK values (see pK_tol, Table S3) Conc2Beta.m Convert a buffer concentration to a BC value PlotBCfit.m Plotting function to generate BC curve plots AdjpKa_AB.m Adjust a pK value for temperature and ionic strength GetAdjC.m Calculate the salt of an acid or base needed to optimize pH (M/L) fmincon.m Matlab function (Optimization toolbox) fminsearch.m Matlab function (Optimization toolbox)

Table 2. Example parameter file

 Parameter Value Description minDeltapH 0.02 Smallest pH change allowed due to a acid or base addition Order 15 Number of terms in the Trigonometric regression MinConc 0.001 Minimim buffer concentration from model predicitons (M/L) WaterSalt 2 NaCl percent for calculating the water BC curve NaClpercent 2 NaCl percent for calculating the BC curve NpKs 7 Initial number of pK values for the BC model fit UB 12 Upper bound for the pH range for BC curve fit LB 2 Lower bound for the pH range for BC curve fit pK_tol 0.2 Minimum pH difference in adjacent pK values trim_beg 0 Number of data points to trim from the start of the BC curve trim_end 10 Number of data points to trim from the end of the BC curve adjC 0 Concentration (M/L) for salt of an acid or base (cations = -val) minGap 0.3 Minimum gap in pH for linear autofilling BC curve increment 0.1 Increment of pH value for gap filling of BC curve