Download 
Welcome to the software download area! Please contact us if you have problems or questions.
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 C_{a} and C_{b} for the acids and bases and the corresponding equilibrium constants K_{a} or K_{b} 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 K_{w} 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^{+}] = (C_{a}K_{a}/(K_{a} + [H^{+}]))  (C_{b}[H^{+}]/([H^{+}] + K_{b})) + K_{w}/[H^{+}]
To solve Equation 1 for [H^{+}], the solution of a 4^{th} 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 A_{ni}, (a positive value for anions) or C_{ti} (negative value for cations) results in:
2) 0 = ∑(C_{ai}K_{ak}/(K_{ai} + [H^{+}]))  ∑(C_{bi}[H^{+}]/([H^{+}] + K_{bi})) + K_{w}/[H^{+}]  [H^{+}] + ∑ (A_{ni} + C_{ti})
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 A_{n }or C_{t} 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 log_{10} 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 “LiveScript” 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 stepwise 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) = B_{o} + A_{1}sin(αx) + B_{1}cos(αx) + A_{2}sin(2αx) + B_{2}cos(2αx) + A_{3}sin(3αx) ...
where F(x) represented the BC value for a given pH (represented as the variable x), A_{i} and B_{i} 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 (P_{i}) of Equation 5 was generated:
6) ∂SS(E)/∂(P_{i}) = 0
An optimized parameter vector (P_{o }= B_{0}, A_{1}, B_{1}, A_{2}, B_{2}, …) for the terms in Equation 4 was then obtained using Matlab’s ‘mldivide’ operator to solve the matrix equation M*P_{o} = Y_{i}, where the vector Y_{i} consisted of the corresponding Y(x) term for each partial derivative equation (implemented as SCBCfit.m, Table 1). Typically, 15 A_{i} and B_{i} parameters (defined as the parameter “Order”, Table 2) were used for the solution. The resulting trigonometric series model with optimized parameters (P_{o}) 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*(∑(C_{i}K_{i}[H^{+}]/([H^{+}] + K_{i})^{2}) + K_{w}/[H^{+}] + [H^{+}])
where C_{i} and K_{i} are the concentration and dissociation constants for weak acids or bases that contribute to buffering in the solution that was originally titrated and K_{w} is the equilibrium constant for water.
A matrix of C_{i} and K_{i} 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 K_{i} 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 concentrationpK 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 (C_{i} 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 A_{n} and C_{t} 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 A_{n} and C_{t }(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. 94201). New York, NY: John Wiley and Sons.
Eubank, R.L. & Speckman P. (1990). Curve fitting by polynomialtrigonomietric regression. Biometrika 77(1):19.
Newbery, A. C. R. (1970). Trigonometric interpolation and curve fitting. Mathematics of Computation, 24(112), 869876.
Longtin M, Price RE, Mishra R, Breidt F. 2020. Modeling the buffer capacity of ingredients in salad dressing products. J Food Sci 85(4):910917.
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):918925.
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 


Download Form
Please complete the form below to track the use of the service we are providing and to help improve our product.OMB 0518  0032 (02/2022)
According to the Paperwork Reduction Act of 1995, an agency cannot conduct or sponsor, and a person is not required to respond to, a collection of information unless it displays a valid OMB control number. The valid OMB control number for this information collection is 05180032. The time required to complete this information is estimated to vary from one to five minutes with an average of three minutes per response, including time for reviewing instructions, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this collection of information, including suggestions for reducing this burden, to Department of Agriculture, Clearance Officer, OIRM, Room 404W, Washington, DC 20250, and to the Office of Information and Regulatory Affairs, Office of Management and Budget.