#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ What a Buffer Is Actually Resisting Science Journaling Club, Volume 1, Issue 3, Spring 2025 Theme: "Equilibrium Chemistry by Computer" ================================================================================ THE QUESTION ================================================================================ Buffer capacity is usually introduced as a formula and left there. Computed properly across the full range of conditions: (a) where is a buffer strongest, (b) how fast does it fail outside its useful window, (c) and how does blood manage to buffer well at a pH 1.3 units away from the dissociation constant that is supposed to set its optimum? ================================================================================ THE MODEL ================================================================================ Buffer capacity, following Van Slyke's 1922 definition, is beta = d(Cb) / d(pH) the number of moles of strong base that must be added to one litre to raise the pH by one unit. We never use the textbook approximation. Instead we write the exact proton condition for the solution and differentiate it. For a solution containing C mol/L of a monoprotic acid HA with dissociation constant Ka, to which Cb mol/L of strong base (NaOH) has been added, charge balance reads [Na+] + [H+] = [OH-] + [A-] Cb + h = Kw/h + C*Ka/(Ka + h) so that Cb(h) = Kw/h - h + C*Ka/(Ka + h) (1) Because pH = -log10(h), we have dh/d(pH) = -ln(10) * h, and therefore beta(h) = ln(10) * [ Kw/h + h + C*Ka*h/(Ka + h)^2 ] (2) Equation (2) is the analytic target. Everything the program reports is checked against it by two independent numerical routes: beta_num a five-point finite difference of Cb(pH) from equation (1); beta_solve add a measured dose of strong base, re-solve the full nonlinear charge balance by bisection for the new pH, and divide. The second route uses no derivative at all. It is the computational equivalent of picking up a burette. For several buffers in the same beaker the buffer terms add, because each conjugate pair takes up protons independently: beta = ln(10) * [ Kw/h + h + SUM_i C_i*Ka_i*h/(Ka_i + h)^2 ] (3) For the carbon dioxide system we treat two cases. CLOSED. Total dissolved carbonate CT is fixed: the beaker is stoppered. With D = h^2 + K1*h + K1*K2, Cb(h) = Kw/h - h + CT*(K1*h + 2*K1*K2)/D (4) beta = ln(10)*[ Kw/h + h + CT*(K1*h^3 + 4*K1*K2*h^2 + K1^2*K2*h)/D^2 ] (5) OPEN. Dissolved carbon dioxide S = alpha * pCO2 is held constant, because the lungs replace or remove it faster than the chemistry drifts. Then [HCO3-] = K1*S/h, [CO3--] = K1*K2*S/h^2 Cb(h) = Kw/h - h + K1*S/h + 2*K1*K2*S/h^2 (6) beta = ln(10)*[ Kw/h + h + K1*S/h + 4*K1*K2*S/h^2 ] (7) Equation (7) has no maximum anywhere in the physiological range. That is the whole answer to question (c), and the program confirms it numerically. ================================================================================ ASSUMPTIONS ================================================================================ 1. Ideal solution. Activities are replaced by concentrations throughout. For the carbon dioxide system we compensate by using APPARENT (concentration) constants measured at plasma ionic strength, which is the standard clinical practice, rather than thermodynamic constants. 2. Equilibrium is instantaneous. No kinetics. The hydration of CO2 to H2CO3 is slow enough to matter in vivo (carbonic anhydrase exists precisely because of it) and we ignore that entirely. 3. [H2CO3] and dissolved CO2 are lumped into one species H2CO3*, as is universal in this field. K1 is therefore an apparent constant. 4. Temperature is fixed: 25 C for the generic buffer work (pKw 14.00), 37 C for the blood work (pKw 13.62). 5. Strong base is added without changing the volume. 6. Dilution on titration is neglected; Cb is in moles per litre of the original solution. 7. In the open system, pCO2 is a fixed boundary condition. Respiration is modelled as an infinitely fast, infinitely capable controller. ================================================================================ LIMITATIONS, STATED PLAINLY ================================================================================ * No solution was ever mixed. The club has no laboratory. Every number below is arithmetic performed on a model of a solution. The computation is the experiment, and where the model is wrong the numbers are wrong with it. * Real blood is not the two-buffer cartoon in section 8. Haemoglobin, plasma protein, phosphate and bone all contribute, and haemoglobin's buffer value changes with oxygenation (the Haldane effect). We model the carbon dioxide system alone and say so wherever we quote a number. * Holding pCO2 perfectly constant is the strongest assumption in the study. A real acid load raises ventilation but does not pin pCO2; a respiratory response that only partly compensates gives a capacity between our open and closed answers. Section 10 quantifies that with a partial-compensation parameter. * Apparent pK2' for plasma is genuinely uncertain, and the reported open/closed ratio moves by several percent across the published range. We carry that uncertainty through a Monte Carlo rather than hiding it. * Activity corrections at ionic strength 0.16 are folded into the apparent constants rather than computed, so the model cannot be pushed to seawater or to concentrated solutions without redoing the constants. ================================================================================ RUN ================================================================================ python buffer-capacity.py > buffer-capacity-output.txt Optional, used to build the figures in the article: python buffer-capacity.py --svg Master seed 20250321. Runtime about 20 seconds on a school laptop. """ import math import os import sys import platform LN10 = math.log(10.0) MASTER_SEED = 20250321 try: import numpy as np HAVE_NUMPY = True NUMPY_VERSION = np.__version__ except Exception: # pragma: no cover np = None HAVE_NUMPY = False NUMPY_VERSION = "not available" # ============================================================================== # CONSTANTS # ============================================================================== PKW_25 = 14.000 # water autoprotolysis, 25 C, conventional value PKW_37 = 13.620 # water autoprotolysis, 37 C KW_25 = 10.0 ** (-PKW_25) KW_37 = 10.0 ** (-PKW_37) PKA_ACETIC = 4.756 # acetic acid, 25 C (Harned & Ehlers 1932) PKA_PHOS2 = 7.198 # H2PO4- / HPO4-- , 25 C PKA_CIT2 = 4.761 # citric acid second step, 25 C PKA_CIT3 = 6.396 # citric acid third step, 25 C PKA_TRIS = 8.072 # tris, 25 C # Blood / plasma, 37 C, apparent (concentration) constants at ionic strength 0.16 PK1_PLASMA = 6.100 # apparent pK1' used in the clinical Henderson-Hasselbalch PK2_PLASMA = 9.800 # apparent pK2' for plasma; genuinely uncertain, see MC ALPHA_CO2 = 0.0307 # CO2 solubility in plasma, mmol/L per mmHg, 37 C PCO2_ART = 40.0 # arterial pCO2, mmHg PH_BLOOD = 7.400 # arterial pH # Published comparison values, quoted in the output for the club's own check. LIT_OPEN_BETA = 2.302585092994046 * 24.0 # mmol/L per pH. The standard teaching # value for the open CO2 system at pH 7.4 is ln(10) times a # round bicarbonate of 24 mmol/L, i.e. 55.26. We quote the # arithmetic rather than a rounded number from a textbook. LIT_PLASMA_NONBICARB = 7.7 # mmol/L per pH, non-bicarbonate buffer value of plasma LIT_BLOOD_NONBICARB = 30.0 # mmol/L per pH, non-bicarbonate buffer value of whole blood # ============================================================================== # CORE CHEMISTRY # ============================================================================== def cb_multi(h, buffers, kw): """Strong base added, mol/L, that puts the solution at [H+] = h. buffers is a list of (total concentration, Ka) pairs. Equation (1)/(3). Negative values mean strong acid was added instead.""" s = kw / h - h for c, ka in buffers: s += c * ka / (ka + h) return s def beta_multi(h, buffers, kw): """Analytic buffer capacity, mol/L per pH unit. Equation (2)/(3).""" s = kw / h + h for c, ka in buffers: s += c * ka * h / (ka + h) ** 2 return LN10 * s def dbeta_dph(h, buffers, kw): """Analytic d(beta)/d(pH). Used to locate the maximum to machine precision.""" s = -kw / h + h for c, ka in buffers: s += c * ka * h * (ka - h) / (ka + h) ** 3 return -LN10 * LN10 * s def cb_carb_closed(h, ct, k1, k2, kw): """Equation (4): stoppered carbonate system, total carbonate fixed.""" d = h * h + k1 * h + k1 * k2 return kw / h - h + ct * (k1 * h + 2.0 * k1 * k2) / d def beta_carb_closed(h, ct, k1, k2, kw): """Equation (5).""" d = h * h + k1 * h + k1 * k2 num = k1 * h ** 3 + 4.0 * k1 * k2 * h * h + k1 * k1 * k2 * h return LN10 * (kw / h + h + ct * num / (d * d)) def cb_carb_open(h, s_co2, k1, k2, kw): """Equation (6): dissolved CO2 pinned by respiration.""" return kw / h - h + k1 * s_co2 / h + 2.0 * k1 * k2 * s_co2 / (h * h) def beta_carb_open(h, s_co2, k1, k2, kw): """Equation (7).""" return LN10 * (kw / h + h + k1 * s_co2 / h + 4.0 * k1 * k2 * s_co2 / (h * h)) def speciation_carb(h, s_co2, k1, k2): """From dissolved CO2 and h, the other two carbonate species.""" hco3 = k1 * s_co2 / h co3 = k1 * k2 * s_co2 / (h * h) return s_co2, hco3, co3 # ============================================================================== # SOLVERS # ============================================================================== def ph_from_base(target_cb, cb_fn, lo=-4.0, hi=19.0, iters=90): """Invert Cb(pH) by bisection. Cb is strictly increasing in pH, so this is unconditionally safe given a bracket. 90 halvings of a 23-unit bracket take the interval far below double precision, so the answer is exact to the last representable bit of pH.""" for _ in range(iters): mid = 0.5 * (lo + hi) if cb_fn(10.0 ** (-mid)) < target_cb: lo = mid else: hi = mid return 0.5 * (lo + hi) def bisect_root(f, lo, hi, iters=90): """Plain bisection for a function that changes sign on [lo, hi].""" flo = f(lo) for _ in range(iters): mid = 0.5 * (lo + hi) fm = f(mid) if (fm < 0.0) == (flo < 0.0): lo = mid flo = fm else: hi = mid return 0.5 * (lo + hi) def beta_five_point(ph, cb_fn, delta=1.0e-3): """Five-point central difference of Cb with respect to pH. Truncation error is order delta^4, i.e. about 1e-12 here.""" def c(p): return cb_fn(10.0 ** (-p)) return (c(ph - 2 * delta) - 8 * c(ph - delta) + 8 * c(ph + delta) - c(ph + 2 * delta)) / (12.0 * delta) def beta_by_titration(ph, cb_fn, target_dph=1.0e-4): """Add and remove a measured dose of strong base, re-solve for pH each time, and divide. No derivative is taken anywhere in this routine. Two passes. The first uses a deliberately tiny dose to get a rough idea of the capacity; the second picks a dose that moves the pH by about target_dph, which keeps the truncation error and the bisection noise in balance. Nothing analytic is consulted at any point.""" cb0 = cb_fn(10.0 ** (-ph)) probe = 1.0e-9 rough = 2.0 * probe / (ph_from_base(cb0 + probe, cb_fn) - ph_from_base(cb0 - probe, cb_fn)) dose = max(rough * target_dph, 1.0e-15) ph_up = ph_from_base(cb0 + dose, cb_fn) ph_dn = ph_from_base(cb0 - dose, cb_fn) return 2.0 * dose / (ph_up - ph_dn) def find_beta_max(buffers, kw, guess, half=0.45): """Locate the capacity maximum by bisecting the analytic derivative inside a bracket that is known to straddle a maximum rather than a minimum. beta has minima as well as maxima once the water walls are included, so a blind bisection over the whole pH range is not safe. We require d(beta)/d(pH) > 0 on the left of the bracket and < 0 on the right, which is a maximum and nothing else. Returns None when no such bracket exists, which happens when the buffer is too weak or too far from neutrality to produce a peak at all.""" def f(p): return dbeta_dph(10.0 ** (-p), buffers, kw) for _ in range(45): lo, hi = guess - half, guess + half if f(lo) > 0.0 > f(hi): return bisect_root(f, lo, hi) half *= 0.5 return None def half_capacity_window(buffers, kw, ph_peak, beta_peak, span=7.0): """pH values either side of the peak where capacity has fallen to half.""" target = 0.5 * beta_peak def f(p): return beta_multi(10.0 ** (-p), buffers, kw) - target lo_edge = ph_peak - span hi_edge = ph_peak + span # The water walls climb again at the extremes, so walk inwards until the # sign really is negative before bisecting. while f(lo_edge) > 0.0 and lo_edge < ph_peak - 0.05: lo_edge += 0.05 while f(hi_edge) > 0.0 and hi_edge > ph_peak + 0.05: hi_edge -= 0.05 lo = bisect_root(f, lo_edge, ph_peak) hi = bisect_root(f, hi_edge, ph_peak) return lo, hi # ============================================================================== # OUTPUT HELPERS # ============================================================================== def rule(ch="="): print(ch * 78) def head(n, title): print() rule() print("%d. %s" % (n, title)) rule() def row(label, value, width=52): print("%-*s: %s" % (width, label, value)) # ============================================================================== # FIGURE DATA, collected as the program runs # ============================================================================== FIG = {} # ============================================================================== # SECTION 0. PROVENANCE AND SOLVER SANITY # ============================================================================== def section_header(): rule() print("WHAT A BUFFER IS ACTUALLY RESISTING") print("Science Journaling Club, computational study") print("Volume 1, Issue 3, Spring 2025, 'Equilibrium Chemistry by Computer'") rule() row("master seed", MASTER_SEED) row("python", platform.python_version()) row("numpy", NUMPY_VERSION) row("rng", "numpy PCG64 seeded from the master seed" if HAVE_NUMPY else "python Mersenne Twister seeded from the master seed") row("pKw at 25 C / 37 C", "%.3f / %.3f" % (PKW_25, PKW_37)) row("units", "concentrations mol/L, beta mol/L per pH unit") print() print("No solution was mixed. The club has no laboratory. The computation") print("is the experiment, and every number below is arithmetic on a model.") def section_sanity(): head(0, "DOES THE SOLVER GET SIMPLE, KNOWN ANSWERS RIGHT?") c = 0.100 ka = 10.0 ** (-PKA_ACETIC) buf = [(c, ka)] def cb(h): return cb_multi(h, buf, KW_25) print("Test solution: %.3f M acetic acid, pKa %.3f, 25 C." % (c, PKA_ACETIC)) print() ph0 = ph_from_base(0.0, cb) short0 = -math.log10(math.sqrt(ka * c)) row("pH of the unbuffered acid, exact solver", "%.6f" % ph0) row("textbook sqrt(Ka*C) shortcut", "%.6f" % short0) row("difference", "%+.6f" % (ph0 - short0)) print() ph_half = ph_from_base(0.5 * c, cb) row("pH at half neutralisation, exact solver", "%.6f" % ph_half) row("pKa", "%.6f" % PKA_ACETIC) row("difference", "%+.6f" % (ph_half - PKA_ACETIC)) print() ph_eq = ph_from_base(c, cb) kb = KW_25 / ka # No dilution in this model, so the conjugate base sits at the full C. short_eq = PKW_25 + math.log10(math.sqrt(kb * c)) row("pH at the equivalence point, exact solver", "%.6f" % ph_eq) row("textbook weak-base shortcut", "%.6f" % short_eq) row("difference", "%+.6f" % (ph_eq - short_eq)) print() worst = 0.0 for k in range(0, 141): p = 0.5 + 0.1 * k back = ph_from_base(cb(10.0 ** (-p)), cb) worst = max(worst, abs(back - p)) row("solver round trip, max |pH_in - pH_out| over 0.5-14.5", "%.3e" % worst) print() print("The shortcuts are not wrong so much as approximate. The solver keeps") print("the terms the shortcuts throw away, which is why it can be trusted at") print("the extremes where the shortcuts fall apart.") # ============================================================================== # SECTION 1. VALIDATION OF THE NUMERICAL DERIVATIVE # ============================================================================== def section_validate_derivative(): head(1, "VALIDATION: NUMERICAL DERIVATIVE AGAINST THE ANALYTIC EXPRESSION") c = 0.100 ka = 10.0 ** (-PKA_ACETIC) buf = [(c, ka)] def cb(h): return cb_multi(h, buf, KW_25) print("System: %.3f M acetic acid, pKa %.3f, 25 C, ideal solution." % (c, PKA_ACETIC)) print("Three independent routes to the same quantity:") print(" beta_exact ln10 [ Kw/h + h + C Ka h/(Ka+h)^2 ] equation (2)") print(" beta_num five-point finite difference of Cb(pH)") print(" beta_solve dose with strong base, re-solve for pH, divide") print() print("%6s %16s %16s %14s %12s %16s %12s" % ("pH", "beta_exact", "beta_num", "diff", "rel", "beta_solve", "rel")) print("-" * 78) max_rel_num = 0.0 max_rel_sol = 0.0 rows = [] p = 1.0 while p <= 13.0001: h = 10.0 ** (-p) be = beta_multi(h, buf, KW_25) bn = beta_five_point(p, cb) bs = beta_by_titration(p, cb) dn = bn - be rn = dn / be rs = (bs - be) / be max_rel_num = max(max_rel_num, abs(rn)) max_rel_sol = max(max_rel_sol, abs(rs)) rows.append((p, be, bn, bs, rn, rs)) print("%6.2f %16.10f %16.10f %+14.3e %+12.2e %16.10f %+12.2e" % (p, be, bn, dn, rn, bs, rs)) p += 0.5 print("-" * 78) row("largest relative difference, finite difference", "%.3e" % max_rel_num) row("largest relative difference, titration route", "%.3e" % max_rel_sol) print() print("VERDICT: all three routes agree. The finite difference matches the") print("analytic expression to the limit of double precision. The titration") print("route, which never differentiates anything, agrees to 1e-7, which is") print("the bisection tolerance on the two re-solved pH values.") FIG["validation"] = rows # ============================================================================== # SECTION 2. WHERE THE MAXIMUM IS # ============================================================================== def section_maximum(): head(2, "VALIDATION: IS THE MAXIMUM REALLY AT pH = pKa?") print("Claim under test: buffer capacity peaks where pH equals pKa, with a") print("peak value of ln(10)*C/4.") print() print("Setting d(beta)/d(pH) = 0 in equation (2) gives") print(" -Kw/h + h + C Ka h (Ka - h)/(Ka + h)^3 = 0") print("At h = Ka the third term vanishes but the first two do not, unless") print("Ka^2 = Kw. So the peak sits at pH = pKa only when pKa = pKw/2.") print("Expanding about h = Ka(1 + eps) gives the displacement") print(" delta_pH ~ -(8/(C ln10)) * (Ka - Kw/Ka)") print() print("%7s %8s %14s %14s %14s %12s" % ("pKa", "C", "pH_max", "pH_max - pKa", "predicted", "found/pred")) print("-" * 78) rows = [] n_peak = 0 n_none = 0 worst_small = 0.0 for pka in (2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0): for c in (0.100, 0.010): ka = 10.0 ** (-pka) buf = [(c, ka)] pmax = find_beta_max(buf, KW_25, pka) pred = -(8.0 / (c * LN10)) * (ka - KW_25 / ka) if pmax is None: n_none += 1 rows.append((pka, c, None, None, pred, None)) print("%7.1f %8.3f %14s %14s %+14.3e %12s" % (pka, c, "no peak", "-", pred, "-")) continue n_peak += 1 disp = pmax - pka ratio = disp / pred if pred != 0.0 else float("nan") rows.append((pka, c, pmax, disp, pred, ratio)) if 4.0 <= pka <= 10.0 and c >= 0.100: worst_small = max(worst_small, abs(disp)) print("%7.1f %8.3f %14.8f %+14.3e %+14.3e %12.5f" % (pka, c, pmax, disp, pred, ratio)) print("-" * 78) row("cases with a genuine local maximum", n_peak) row("cases with no local maximum at all", n_none) row("largest |pH_max - pKa| for pKa 4-10 at 0.100 M", "%.3e pH units" % worst_small) print() print("VERDICT: the textbook claim holds to within 0.0036 pH units for any") print("0.100 M buffer between pKa 4 and pKa 10, and to within 0.00035 pH") print("units between pKa 5 and pKa 9.") print("The residual is not noise. It is the water term, and the first-order") print("prediction reproduces it with a found/predicted ratio of 1.000 wherever") print("the displacement is small enough for the expansion to hold. At") print("pKa 7.000 the displacement is identically zero, because that is the one") print("place where Ka = Kw/Ka. The club did not expect a displacement at all") print("and spent an afternoon looking for a bug before finding the algebra.") print() print("The rows marked 'no peak' are the real surprise. A pKa 2 or pKa 12") print("buffer at these concentrations produces no maximum whatsoever: the") print("capacity falls monotonically away from the strong-acid wall, and the") print("buffer contributes a shoulder rather than a peak. Below about") print("0.01 M the same thing happens at pKa 3. Anything you buy labelled as") print("a buffer for pH 2 is relying on the solvent, not on the conjugate pair.") print() print("Peak height against the analytic prediction ln(10)*C/4:") print() print("%7s %8s %16s %16s %14s %12s" % ("pKa", "C", "beta_max", "ln10*C/4", "difference", "rel")) print("-" * 78) peak_rows = [] for pka in (3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0): for c in (0.100, 0.010): ka = 10.0 ** (-pka) buf = [(c, ka)] pmax = find_beta_max(buf, KW_25, pka) if pmax is None: continue bmax = beta_multi(10.0 ** (-pmax), buf, KW_25) pred = LN10 * c / 4.0 water = LN10 * (KW_25 / 10.0 ** (-pmax) + 10.0 ** (-pmax)) peak_rows.append((pka, c, pmax, bmax, pred, water)) print("%7.1f %8.3f %16.10f %16.10f %+14.3e %12.6f" % (pka, c, bmax, pred, bmax - pred, (bmax - pred) / pred)) print("-" * 78) print() print("VERDICT: ln(10)*C/4 is exact for the buffer term and low for the total") print("by very close to the water contribution ln(10)*(Kw/h + h) at the peak.") print("For pKa 7 at 0.100 M that correction is eight parts per million. For") print("pKa 3 at 0.100 M it is 4%, because there the solvent is helping.") return peak_rows # ============================================================================== # SECTION 3. HOW FAST IT FAILS # ============================================================================== def section_window(): head(3, "HOW FAST CAPACITY FALLS OUTSIDE THE WINDOW") print("Ignoring water, beta/beta_max = 4r/(1+r)^2 with r = 10^(pKa - pH).") print("Setting that to one half gives r^2 - 6r + 1 = 0, so r = 3 + 2*sqrt(2)") print("and the half-capacity points sit at pKa +/- log10(3 + 2*sqrt(2)).") print() analytic_half = math.log10(3.0 + 2.0 * math.sqrt(2.0)) row("analytic half-width, pH units", "%.9f" % analytic_half) row("analytic full width", "%.9f" % (2.0 * analytic_half)) print() print("Retained fraction of peak capacity at a given distance from pKa,") print("water excluded, 4r/(1+r)^2 with r = 10^d:") print() print("%10s %16s %14s" % ("d (pH)", "fraction", "1 in")) print("-" * 78) retained = [] for d in (0.0, 0.25, 0.5, 0.7655, 1.0, 1.5, 2.0, 2.5, 3.0, 4.0): r = 10.0 ** d f = 4.0 * r / (1.0 + r) ** 2 retained.append((d, f)) print("%10.4f %16.9f %14.1f" % (d, f, 1.0 / f)) print("-" * 78) print() print("Now the same thing computed, not derived, from the full model with") print("water present. 0.100 M buffers, 25 C:") print() print("%7s %12s %12s %12s %12s %14s" % ("pKa", "pH_max", "half lo", "half hi", "width", "width - exact")) print("-" * 78) win_rows = [] for pka in (4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0): buf = [(0.100, 10.0 ** (-pka))] pmax = find_beta_max(buf, KW_25, pka) bmax = beta_multi(10.0 ** (-pmax), buf, KW_25) lo, hi = half_capacity_window(buf, KW_25, pmax, bmax) w = hi - lo win_rows.append((pka, pmax, lo, hi, w)) print("%7.1f %12.6f %12.6f %12.6f %12.6f %+14.3e" % (pka, pmax, lo, hi, w, w - 2.0 * analytic_half)) print("-" * 78) print() print("VERDICT: the computed half-capacity window matches the analytic") print("1.531103 pH units to five parts in a hundred thousand at pKa 7, to") print("two parts in ten thousand at pKa 6 and 8, and to 1.8% at pKa 4 and") print("10, where the water walls have begun to crowd the window. Nothing") print("here is noise; every deviation is the solvent. The conventional") print("'pKa plus or minus one' is generous: at one unit out you are down to") print("33.06% of peak, not 50%.") FIG["retained"] = retained return analytic_half, win_rows, retained # ============================================================================== # SECTION 4. CONCENTRATION SCALING # ============================================================================== def section_concentration(): head(4, "DOES CAPACITY SCALE WITH CONCENTRATION?") print("Buffer at pKa 7.000, 25 C. The buffer term of equation (2) is linear") print("in C, so beta_max should be too, once the water floor is subtracted.") print() print("%12s %16s %16s %16s %12s" % ("C (mol/L)", "beta_max", "beta_max/C", "water floor", "buffer/total")) print("-" * 78) rows = [] for c in (0.0001, 0.001, 0.005, 0.010, 0.050, 0.100, 0.250, 0.500, 1.000): buf = [(c, 1.0e-7)] pmax = find_beta_max(buf, KW_25, 7.0) h = 10.0 ** (-pmax) bmax = beta_multi(h, buf, KW_25) water = LN10 * (KW_25 / h + h) rows.append((c, bmax, water)) print("%12.4f %16.10f %16.9f %16.10f %12.6f" % (c, bmax, bmax / c, water, (bmax - water) / bmax)) print("-" * 78) print() row("ln(10)/4, the predicted slope", "%.9f" % (LN10 / 4.0)) print() print("VERDICT: beta_max/C converges on ln(10)/4 = 0.575646 from above, the") print("excess being the water floor of 4.605e-7 mol/L per pH unit at pH 7.") print("At 0.1 mM the water contributes 0.8% of the total. At 1 M it") print("contributes eight parts in ten million. A very dilute buffer is not") print("really a buffer; it is water with a flavour.") return rows # ============================================================================== # SECTION 5. TWO COMPONENTS # ============================================================================== def section_two_component(): head(5, "TWO BUFFERS IN THE SAME BEAKER") print("Total concentration fixed at 0.100 M, split evenly between two") print("monoprotic buffers whose pKa values straddle pH 7.000 by +/- D/2.") print() print("The pair produces one flat-topped peak while D is small and two") print("separate peaks once D is large. The changeover is where the second") print("derivative at the midpoint vanishes, and that condition reduces to") print("cosh(ln10 * D/2) = 2, so") print(" D* = 2*arccosh(2)/ln(10)") d_star = 2.0 * math.acosh(2.0) / LN10 row("analytic split threshold D*", "%.9f" % d_star) print() print("Width below is the span between the outermost pH values where the") print("mixture still delivers half of what a single 0.100 M buffer manages") print("at its own peak, i.e. beta >= ln(10)*0.100/8 = %.8f." % (LN10 * 0.1 / 8.0)) print() thr = LN10 * 0.100 / 8.0 print("%8s %16s %16s %14s %8s %12s" % ("D", "beta(pH 7)", "beta_max", "peak at", "shape", "width")) print("-" * 78) rows = [] for d in (0.0, 0.5, 0.9, 1.0, d_star, 1.2, 1.5, 2.0, 2.5, 3.0): pk1 = 7.0 - d / 2.0 pk2 = 7.0 + d / 2.0 buf = [(0.050, 10.0 ** (-pk1)), (0.050, 10.0 ** (-pk2))] bmid = beta_multi(1.0e-7, buf, KW_25) pk = find_beta_max(buf, KW_25, pk1, half=0.45) if pk is None or d < d_star: pk = find_beta_max(buf, KW_25, 7.0, half=0.45) best_b = beta_multi(10.0 ** (-pk), buf, KW_25) if pk is not None else bmid best_p = pk if pk is not None else 7.0 shape = "single" if abs(best_p - 7.0) < 1e-6 else "split" def f(p, bb=buf): return beta_multi(10.0 ** (-p), bb, KW_25) - thr n = 3001 grid = [3.0 + 8.0 * i / (n - 1) for i in range(n)] above = [p for p in grid if f(p) > 0.0] if above: lo = bisect_root(f, above[0] - 8.0 / (n - 1), above[0]) hi = bisect_root(f, above[-1] + 8.0 / (n - 1), above[-1]) width = hi - lo else: width = 0.0 rows.append((d, bmid, best_b, best_p, shape, width)) print("%8.4f %16.10f %16.10f %14.6f %8s %12.6f" % (d, bmid, best_b, best_p, shape, width)) print("-" * 78) print() print("VERDICT: the numerical scan flips from a single central peak to a") print("split pair between D = %.4f and D = 1.2000, which brackets the" % d_star) print("analytic threshold %.6f. Spacing the pKa values by roughly that" % d_star) print("amount buys the widest useful plateau. A single 0.100 M buffer holds") print("half its own peak over 1.531 pH units. The best-spaced pair holds") print("half of that same reference capacity over about two and a third pH") print("units, for a peak that is lower. Width is bought with height.") print() print("Two real pairs, each at 0.050 M per component, 25 C:") print() pairs = [ ("acetate + phosphate", PKA_ACETIC, PKA_PHOS2), ("citrate 2nd + 3rd", PKA_CIT2, PKA_CIT3), ("phosphate + tris", PKA_PHOS2, PKA_TRIS), ] real_rows = [] for name, a, b in pairs: buf = [(0.050, 10.0 ** (-a)), (0.050, 10.0 ** (-b))] sep = b - a # find the deepest dip between the two pKa values dip_p, dip_b = None, float("inf") n = 2001 for i in range(n): p = a + (b - a) * i / (n - 1) bb = beta_multi(10.0 ** (-p), buf, KW_25) if bb < dip_b: dip_b, dip_p = bb, p top_p, top_b = None, -1.0 for i in range(3001): p = 2.0 + 10.0 * i / 3000.0 bb = beta_multi(10.0 ** (-p), buf, KW_25) if bb > top_b: top_b, top_p = bb, p split = (b - a) > d_star real_rows.append((name, a, b, sep, dip_p, dip_b, top_b, split)) print(" %-22s pKa %.3f and %.3f, separation %.3f" % (name, a, b, sep)) print(" highest capacity anywhere : %.6f at pH %.3f" % (top_b, top_p)) print(" lowest point between them : %.6f at pH %.3f" % (dip_b, dip_p)) print(" that dip, relative to peak: %.4f" % (dip_b / top_b)) print(" past the split threshold : %s" % ("yes" if split else "no")) print() print("Acetate and phosphate sit 2.442 pH units apart, well past the") print("threshold, so the mixture sags to %.1f%% of its peak in the middle." % (100.0 * real_rows[0][5] / real_rows[0][6])) print("Citrate's own second and third steps sit 1.635 apart and sag only to") print("%.1f%%. Phosphate and tris sit 0.874 apart, inside the threshold, and" % (100.0 * real_rows[1][5] / real_rows[1][6])) print("show no dip in the middle at all. Their capacity is one broad peak") print("at pH 7.637, and the lowest value between the two pKa figures sits at") print("an end of that interval rather than inside it.") return d_star, rows, real_rows # ============================================================================== # SECTION 6. THE CARBON DIOXIDE SYSTEM, CLOSED AND OPEN # ============================================================================== def blood_state(ph=PH_BLOOD, pco2=PCO2_ART, pk1=PK1_PLASMA, pk2=PK2_PLASMA, alpha=ALPHA_CO2, kw=KW_37): """Everything the model knows about a blood-like solution at one setting.""" h = 10.0 ** (-ph) k1 = 10.0 ** (-pk1) k2 = 10.0 ** (-pk2) s = alpha * pco2 / 1000.0 # mmol/L per mmHg -> mol/L co2, hco3, co3 = speciation_carb(h, s, k1, k2) ct = co2 + hco3 + co3 b_open = beta_carb_open(h, s, k1, k2, kw) b_closed = beta_carb_closed(h, ct, k1, k2, kw) return { "h": h, "k1": k1, "k2": k2, "s": s, "co2": co2, "hco3": hco3, "co3": co3, "ct": ct, "open": b_open, "closed": b_closed, "ratio": b_open / b_closed, } def section_carbonate(): head(6, "THE CARBON DIOXIDE SYSTEM AT BLOOD pH, CLOSED AND OPEN") st = blood_state() print("Conditions: 37 C, pH %.3f, pCO2 %.1f mmHg, apparent pK1' %.3f," % (PH_BLOOD, PCO2_ART, PK1_PLASMA)) print("apparent pK2' %.3f, CO2 solubility %.4f mmol/L per mmHg, pKw %.3f." % (PK2_PLASMA, ALPHA_CO2, PKW_37)) print() print("Speciation derived from those inputs, not assumed:") row(" dissolved CO2 (H2CO3*)", "%9.4f mmol/L" % (1000 * st["co2"])) row(" bicarbonate HCO3-", "%9.4f mmol/L" % (1000 * st["hco3"])) row(" carbonate CO3--", "%9.4f mmol/L" % (1000 * st["co3"])) row(" total carbonate CT", "%9.4f mmol/L" % (1000 * st["ct"])) row(" clinical value for HCO3- at these settings", "24 mmol/L") print() print("The model was given pH and pCO2 and produced %.2f mmol/L bicarbonate." % (1000 * st["hco3"])) print("That is the number a blood gas analyser reports for a healthy adult,") print("and nothing in the calculation was tuned to make it come out.") print() row("beta, CLOSED system (CT fixed)", "%12.6f mol/L per pH" % st["closed"]) row(" same, in clinical units", "%12.3f mmol/L per pH" % (1000 * st["closed"])) row("beta, OPEN system (pCO2 fixed)", "%12.6f mol/L per pH" % st["open"]) row(" same, in clinical units", "%12.3f mmol/L per pH" % (1000 * st["open"])) row("RATIO open / closed", "%12.4f" % st["ratio"]) print() approx = 1.0 + 10.0 ** (PH_BLOOD - PK1_PLASMA) row("simple prediction 1 + 10^(pH - pK1')", "%12.4f" % approx) row("difference from the full model", "%+12.4f" % (st["ratio"] - approx)) row("relative difference", "%+12.4f" % ((st["ratio"] - approx) / approx)) print() print("Where the simple prediction comes from: drop the carbonate ion and") print("the water term, and beta_closed = ln10*CT*a0*a1 while beta_open =") print("ln10*CT*a1, so the ratio is 1/a0 = 1 + K1/h. The full model sits") print("6.9% below that because of the carbonate ion, which the shortcut") print("drops. Carbonate contributes four times its own concentration to the") print("closed capacity through the 4*K1*K2*h^2 term and rather less to the") print("open one, so keeping it raises the denominator more than the") print("numerator.") print() print("Published comparison, for the club's own check:") row(" teaching value ln(10) x 24 mmol/L bicarbonate", "%.2f mmol/L per pH" % LIT_OPEN_BETA) row(" this model", "%.1f mmol/L per pH" % (1000 * st["open"])) row(" difference", "%+.2f mmol/L per pH" % (1000 * st["open"] - LIT_OPEN_BETA)) row(" relative", "%+.1f%%" % (100.0 * (1000 * st["open"] - LIT_OPEN_BETA) / LIT_OPEN_BETA)) row(" measured plasma buffer capacity, all buffers", "16.1 mmol/L per pH") row(" measured whole-blood buffer capacity, all buffers", "38.5 mmol/L per pH") print() print("VERDICT: the open system buffers %.1f times better than the closed" % st["ratio"]) print("one at blood pH. That reproduces the standard physiological result,") print("which is usually quoted as roughly twentyfold. The model's open-system") print("capacity sits %.1f%% above the ln(10) x 24 teaching value, the excess being" % (100.0 * (1000 * st["open"] - LIT_OPEN_BETA) / LIT_OPEN_BETA)) print("the half millimole of extra bicarbonate our inputs imply plus the carbonate") print("ion the teaching version drops.") return st # ============================================================================== # SECTION 7. THE OPEN SYSTEM HAS NO OPTIMUM # ============================================================================== def section_no_optimum(st): head(7, "THE OPEN SYSTEM HAS NO OPTIMUM AT ALL") print("A closed buffer has a maximum at pH = pKa and falls away on both") print("sides. Equation (7) does not. At fixed pCO2, bicarbonate rises by a") print("factor of ten for every pH unit, so capacity rises with it.") print() print("%8s %16s %16s %14s %16s" % ("pH", "beta_open", "beta_closed", "ratio", "HCO3- (mmol/L)")) print("-" * 78) k1, k2, s = st["k1"], st["k2"], st["s"] rows = [] p = 6.0 while p <= 8.4001: h = 10.0 ** (-p) co2, hco3, co3 = speciation_carb(h, s, k1, k2) ct = co2 + hco3 + co3 bo = beta_carb_open(h, s, k1, k2, KW_37) bc = beta_carb_closed(h, ct, k1, k2, KW_37) rows.append((p, bo, bc, hco3)) print("%8.2f %16.8f %16.8f %14.4f %16.4f" % (p, bo, bc, bo / bc, 1000 * hco3)) p += 0.2 print("-" * 78) FIG["carb"] = rows print() print("Monotonicity check. beta_open has exactly one turning point, a") print("minimum, where the strong-acid term h stops dominating and the") print("bicarbonate term K1*S/h takes over. Above that point it rises") print("without limit and never turns back.") turn = None prev = None p = 2.0 while p <= 11.0001: b = beta_carb_open(10.0 ** (-p), s, k1, k2, KW_37) if prev is not None and b > prev and turn is None: turn = p prev = b p += 0.001 increasing = True bad = None prev = None p = turn + 0.01 if turn else 5.0 while p <= 11.0001: b = beta_carb_open(10.0 ** (-p), s, k1, k2, KW_37) if prev is not None and b <= prev: increasing = False bad = p prev = b p += 0.001 row("location of the single minimum", "pH %.3f" % turn) row("beta_open strictly increasing above it, up to pH 11", "yes" if increasing else "no, first decrease at pH %.3f" % bad) row("a closed buffer, by contrast", "has a maximum and falls away both sides") print() print("A closed carbonate system of the same total carbonate peaks at its") print("own pK1' and is worth:") ct = st["ct"] buf_like = beta_carb_closed(10.0 ** (-PK1_PLASMA), ct, k1, k2, KW_37) row(" beta_closed at its own optimum pH %.2f" % PK1_PLASMA, "%.3f mmol/L per pH" % (1000 * buf_like)) row(" beta_open at blood pH 7.40, 1.3 units off", "%.3f mmol/L per pH" % (1000 * st["open"])) row(" open, off-optimum, beats closed, on-optimum, by", "%.2f times" % (st["open"] / buf_like)) print() print("This is the answer to the question the study set out with. Blood") print("does not buffer well at pH 7.4 despite the carbonate pK1' being 6.1.") print("It buffers well BECAUSE the system is open, and openness converts the") print("1.3-unit offset from a penalty into an advantage: the further above") print("pK1' you sit, the more bicarbonate a fixed pCO2 generates.") print() conc_equiv = 4.0 * st["open"] / LN10 row("closed buffer at pKa 7.40 needed to match the open system", "%.1f mmol/L" % (1000 * conc_equiv)) row("total carbonate actually present", "%.1f mmol/L" % (1000 * ct)) row("factor", "%.2f" % (conc_equiv / ct)) return buf_like, conc_equiv # ============================================================================== # SECTION 8. WHAT A DOSE OF ACID ACTUALLY DOES # ============================================================================== def section_dose(st): head(8, "WHAT A DOSE OF ACID ACTUALLY DOES") print("Capacity is a derivative, so it only describes small doses exactly.") print("Here we add strong acid to one litre and re-solve the full equations,") print("with no linearisation anywhere.") print() k1, k2, s, ct = st["k1"], st["k2"], st["s"], st["ct"] def cb_open(h): return cb_carb_open(h, s, k1, k2, KW_37) def cb_closed(h): return cb_carb_closed(h, ct, k1, k2, KW_37) cb0_open = cb_open(st["h"]) cb0_closed = cb_closed(st["h"]) print("%12s %14s %14s %14s %14s" % ("acid mmol/L", "pH open", "pH closed", "dpH open", "dpH closed")) print("-" * 78) rows = [] for dose_mm in (1.0, 2.0, 5.0, 10.0, 20.0): d = dose_mm / 1000.0 po = ph_from_base(cb0_open - d, cb_open) pc = ph_from_base(cb0_closed - d, cb_closed) rows.append((dose_mm, po, pc)) print("%12.1f %14.4f %14.4f %+14.4f %+14.4f" % (dose_mm, po, pc, po - PH_BLOOD, pc - PH_BLOOD)) print("-" * 78) print() d = 0.010 po = ph_from_base(cb0_open - d, cb_open) pc = ph_from_base(cb0_closed - d, cb_closed) lin_o = PH_BLOOD - d / st["open"] lin_c = PH_BLOOD - d / st["closed"] print("Linear prediction against the exact solve, 10 mmol/L of strong acid:") row(" open, exact", "%.4f" % po) row(" open, pH - dose/beta", "%.4f" % lin_o) row(" open, error of the linear estimate", "%+.4f pH" % (lin_o - po)) row(" closed, exact", "%.4f" % pc) row(" closed, pH - dose/beta", "%.4f" % lin_c) row(" closed, error of the linear estimate", "%+.4f pH" % (lin_c - pc)) print() print("VERDICT: for the open system the derivative is an honest guide even") print("at 10 mmol/L, because capacity barely changes over that interval.") print("For the closed system the linear estimate is badly wrong, and wrong") print("in the reassuring direction: it predicts a pH the solution never") print("reaches, because capacity collapses as the bicarbonate is consumed.") print("A metabolic acid load of 10 mmol/L is survivable in the open system") print("and catastrophic in the closed one.") return rows, (po, pc, lin_o, lin_c) # ============================================================================== # SECTION 9. MONTE CARLO ON THE RATIO # ============================================================================== def section_monte_carlo(): head(9, "MONTE CARLO: HOW WELL DO WE KNOW THAT RATIO?") n_trials = 200000 print("The open/closed ratio depends on five inputs the club does not know") print("exactly. We draw them %d times and recompute the ratio each time." % n_trials) print("The distributions below are the club's stated priors. They are") print("assumptions, not measurements, and they are the weakest link in the") print("uncertainty estimate.") print() print(" pH normal(7.400, 0.020), clipped to [7.30, 7.50]") print(" pCO2 normal(40.0, 3.0) mmHg, clipped to [30, 55]") print(" pK1' normal(6.100, 0.020)") print(" pK2' uniform(9.80, 10.40), the published range for plasma") print(" alpha normal(0.0307, 0.0008) mmol/L per mmHg") print() row("master seed", MASTER_SEED) row("trials", n_trials) if HAVE_NUMPY: ss = np.random.SeedSequence(MASTER_SEED) rng = np.random.default_rng(ss) row("generator", "numpy %s, PCG64" % NUMPY_VERSION) ph = np.clip(rng.normal(PH_BLOOD, 0.020, n_trials), 7.30, 7.50) pco2 = np.clip(rng.normal(PCO2_ART, 3.0, n_trials), 30.0, 55.0) pk1 = rng.normal(PK1_PLASMA, 0.020, n_trials) pk2 = rng.uniform(9.80, 10.40, n_trials) alpha = rng.normal(ALPHA_CO2, 0.0008, n_trials) h = 10.0 ** (-ph) k1 = 10.0 ** (-pk1) k2 = 10.0 ** (-pk2) s = alpha * pco2 / 1000.0 hco3 = k1 * s / h co3 = k1 * k2 * s / (h * h) ct = s + hco3 + co3 b_open = LN10 * (KW_37 / h + h + hco3 + 4.0 * co3) d = h * h + k1 * h + k1 * k2 num = k1 * h ** 3 + 4.0 * k1 * k2 * h * h + k1 * k1 * k2 * h b_closed = LN10 * (KW_37 / h + h + ct * num / (d * d)) ratio = b_open / b_closed vals = ratio opens = b_open mean = float(vals.mean()) sd = float(vals.std(ddof=1)) omean = float(opens.mean()) osd = float(opens.std(ddof=1)) qs = [float(x) for x in np.quantile(vals, [0.025, 0.25, 0.50, 0.75, 0.975])] running = np.cumsum(vals) / np.arange(1, n_trials + 1) run_sq = np.cumsum(vals * vals) idx = np.unique(np.round(np.logspace(1, math.log10(n_trials), 140)).astype(int)) idx = idx[idx >= 10] conv = [] for i in idx: m = float(running[i - 1]) var = (float(run_sq[i - 1]) - i * m * m) / (i - 1) conv.append((int(i), m, math.sqrt(max(var, 0.0) / i))) else: # pragma: no cover import random rnd = random.Random(MASTER_SEED) row("generator", "python Mersenne Twister (numpy unavailable)") tot = 0.0 tot2 = 0.0 vals_list = [] conv = [] checkpoints = set(int(round(x)) for x in [10 * (n_trials / 10.0) ** (j / 139.0) for j in range(140)]) for i in range(1, n_trials + 1): ph = min(7.50, max(7.30, rnd.gauss(PH_BLOOD, 0.020))) pco2 = min(55.0, max(30.0, rnd.gauss(PCO2_ART, 3.0))) pk1 = rnd.gauss(PK1_PLASMA, 0.020) pk2 = rnd.uniform(9.80, 10.40) alpha = rnd.gauss(ALPHA_CO2, 0.0008) stx = blood_state(ph, pco2, pk1, pk2, alpha, KW_37) v = stx["ratio"] vals_list.append(v) tot += v tot2 += v * v if i in checkpoints and i >= 10: m = tot / i var = (tot2 - i * m * m) / (i - 1) conv.append((i, m, math.sqrt(max(var, 0.0) / i))) vals_list.sort() mean = tot / n_trials sd = math.sqrt((tot2 - n_trials * mean * mean) / (n_trials - 1)) qs = [vals_list[int(q * (n_trials - 1))] for q in (0.025, 0.25, 0.5, 0.75, 0.975)] omean = float("nan") osd = float("nan") se = sd / math.sqrt(n_trials) print() row("mean ratio", "%.4f" % mean) row("standard deviation across draws", "%.4f" % sd) row("standard error of the mean", "%.5f" % se) row("relative standard error", "%.4f %%" % (100.0 * se / mean)) print() print("%14s %12s %12s %12s %12s" % ("2.5%", "25%", "median", "75%", "97.5%")) print("-" * 78) print("%14.4f %12.4f %12.4f %12.4f %12.4f" % tuple(qs)) print("-" * 78) if HAVE_NUMPY: print() row("mean open-system beta", "%.3f mmol/L per pH" % (1000 * omean)) row("sd of open-system beta", "%.3f mmol/L per pH" % (1000 * osd)) row("standard error", "%.4f mmol/L per pH" % (1000 * osd / math.sqrt(n_trials))) print() print("Convergence of the running mean (every tenth checkpoint shown):") print() print("%12s %14s %14s %14s" % ("trials", "running mean", "running SE", "mean +/- 1 SE")) print("-" * 78) for j, (i, m, s_) in enumerate(conv): if j % 10 == 0 or j == len(conv) - 1: print("%12d %14.4f %14.5f %14s" % (i, m, s_, "%.4f-%.4f" % (m - s_, m + s_))) print("-" * 78) FIG["conv"] = conv print() base = blood_state()["ratio"] zed = (mean - base) / se row("central-value ratio from section 6", "%.4f" % base) row("Monte Carlo mean", "%.4f" % mean) row("difference in standard errors", "%.2f" % zed) print() print("The Monte Carlo mean sits %.0f standard errors from the central-value" % abs(zed)) print("calculation of section 6. That is not a disagreement about chemistry.") print("The central value uses pK2' = 9.800, which is the BOTTOM edge of the") print("uniform prior; the prior's own mean is 10.100. Rerunning section 6 at") print("pK2' = 10.100 gives:") row(" central-value ratio at pK2' = 10.100", "%.4f" % blood_state(pk2=10.10)["ratio"]) print("which sits two percent of one standard deviation from the Monte Carlo") print("mean. The remaining sliver is genuine curvature: the ratio is") print("nonlinear in pH and pK1', so the average of the function is not the") print("function of the average. The standard error, %.5f, measures how many" % se) print("draws we took. It says nothing about how well anyone knows plasma") print("chemistry. The standard deviation, %.3f, is the honest width, and it" % sd) print("is four hundred times larger.") return mean, sd, se, qs, conv, zed # ============================================================================== # SECTION 10. SENSITIVITY AND PARTIAL COMPENSATION # ============================================================================== def section_sensitivity(st): head(10, "SENSITIVITY: WHICH CHOICES WOULD HAVE CHANGED THE ANSWER?") base = st["ratio"] print("One input moved at a time, everything else held at the base case.") print() print("%-34s %12s %12s %12s" % ("variant", "ratio", "change", "beta_open")) print("-" * 78) variants = [ ("base case", dict()), ("pK1' 6.03 (Siggaard-Andersen low)", dict(pk1=6.03)), ("pK1' 6.17", dict(pk1=6.17)), ("pK2' 10.10 (mean of the MC prior)", dict(pk2=10.10)), ("pK2' 10.33 (thermodynamic 25 C)", dict(pk2=10.33)), ("pK2' 9.60", dict(pk2=9.60)), ("pCO2 30 mmHg (hyperventilating)", dict(pco2=30.0)), ("pCO2 55 mmHg (retaining)", dict(pco2=55.0)), ("pH 7.25 (acidotic)", dict(ph=7.25)), ("pH 7.55 (alkalotic)", dict(ph=7.55)), ("alpha 0.0301", dict(alpha=0.0301)), ("alpha 0.0313", dict(alpha=0.0313)), ("pKw 14.00 (25 C water)", dict(kw=KW_25)), ] rows = [] for name, kw_args in variants: v = blood_state(**kw_args) rows.append((name, v["ratio"], v["ratio"] - base, 1000 * v["open"])) print("%-34s %12.4f %+12.4f %12.3f" % (name, v["ratio"], v["ratio"] - base, 1000 * v["open"])) print("-" * 78) print() print("VERDICT: pH and pK1' dominate, as the 1 + 10^(pH - pK1') form") print("predicts. Nothing else moves the ratio by more than a few percent.") print("Dropping to 25 C water changes it in the fifth decimal, which is the") print("clearest possible statement that water itself buffers nothing here.") print() print("HOW GOOD DOES THE LUNG HAVE TO BE?") print("A perfect lung pins pCO2. A stoppered beaker lets pCO2 climb as acid") print("is added, because the acid converts bicarbonate straight back into") print("dissolved carbon dioxide. Let g be the number of decades pCO2 rises") print("per unit fall in pH:") print(" pCO2(pH) = pCO2_0 * 10^(-g*(pH - 7.400))") print("g = 0 is the fully open system of section 6. Negative g is active") print("hyperventilation, which blows carbon dioxide off faster than the") print("chemistry makes it. Somewhere near g = 1 the system behaves like the") print("stoppered beaker. The capacity below is a five-point finite") print("difference of the full proton condition with pCO2 varying this way.") print() k1, k2 = st["k1"], st["k2"] alpha = ALPHA_CO2 def cb_partial(h, g): p = -math.log10(h) pco2 = PCO2_ART * 10.0 ** (-g * (p - PH_BLOOD)) s = alpha * pco2 / 1000.0 return KW_37 / h - h + k1 * s / h + 2.0 * k1 * k2 * s / (h * h) def beta_g(g): return beta_five_point(PH_BLOOD, lambda h: cb_partial(h, g), delta=1.0e-3) print("%10s %18s %16s %18s" % ("g", "beta", "mmol/L per pH", "fraction of open")) print("-" * 78) prows = [] for g in (-0.30, -0.10, 0.0, 0.05, 0.10, 0.25, 0.50, 0.75, 0.90, 1.00): b = beta_g(g) prows.append((g, b)) print("%10.2f %18.8f %16.3f %18.4f" % (g, b, 1000 * b, b / st["open"])) print("-" * 78) print() g_closed = bisect_root(lambda g: beta_g(g) - st["closed"], 0.0, 1.2) row("g that reproduces the closed-system capacity", "%.6f" % g_closed) row("beta at that g", "%.3f mmol/L per pH" % (1000 * beta_g(g_closed))) row("closed-system beta from section 6", "%.3f mmol/L per pH" % (1000 * st["closed"])) print() print("VERDICT: capacity falls almost exactly linearly in g, because the") print("bicarbonate term of equation (7) carries a factor (1 - g). Letting") print("pCO2 drift by a tenth of a decade per pH unit costs %.0f%% of the" % (100.0 * (1.0 - prows[4][1] / st["open"]))) print("capacity. Letting it drift by %.3f decades, which is what a sealed" % g_closed) print("container does, costs %.1f%%. The lung does not have to be perfect." % (100.0 * (1.0 - st["closed"] / st["open"]))) print("It has to be better than a cork.") return rows, prows, g_closed # ============================================================================== # SECTION 11. SUMMARY TABLE # ============================================================================== def section_summary(st, analytic_half, d_star, mc): head(11, "SUMMARY OF HEADLINE NUMBERS") mean, sd, se, qs, conv, zed = mc items = [ ("peak capacity of a monoprotic buffer", "ln(10)*C/4", "%.6f mol/L per pH at C = 0.100" % (LN10 * 0.1 / 4.0)), ("position of that peak", "pH = pKa", "within 0.0036 pH for pKa 4-10 at 0.100 M"), ("exact displacement at pKa 7.000", "zero", "%.3e pH units" % abs(find_beta_max([(0.1, 1e-7)], KW_25, 7.0) - 7.0)), ("half-capacity half-width", "log10(3 + 2 sqrt 2)", "%.6f pH units" % analytic_half), ("capacity retained one pH unit from pKa", "40/121", "%.4f" % (40.0 / 121.0)), ("capacity retained two pH units from pKa", "400/10201", "%.5f" % (400.0 / 10201.0)), ("two-buffer split threshold", "2 arccosh(2)/ln10", "%.6f pH units" % d_star), ("blood, bicarbonate from pH and pCO2", "model output", "%.2f mmol/L" % (1000 * st["hco3"])), ("blood, open-system capacity", "model output", "%.2f mmol/L per pH" % (1000 * st["open"])), ("blood, closed-system capacity", "model output", "%.3f mmol/L per pH" % (1000 * st["closed"])), ("open / closed ratio, central value", "model output", "%.3f" % st["ratio"]), ("open / closed ratio, Monte Carlo", "200000 draws", "%.3f +/- %.4f (SE), sd %.3f" % (mean, se, sd)), ("simple prediction for that ratio", "1 + 10^(pH - pK1')", "%.3f" % (1.0 + 10.0 ** (PH_BLOOD - PK1_PLASMA))), ] for a, b, c in items: print(" %-42s %-22s %s" % (a, b, c)) print() rule() print("END OF OUTPUT") rule() # ============================================================================== # SVG FIGURE EMISSION (optional, for the article) # ============================================================================== class Ax(object): def __init__(self, x0, y0, w, h, xlim, ylim, xlog=False, ylog=False): self.x0, self.y0, self.w, self.h = x0, y0, w, h self.xa, self.xb = xlim self.ya, self.yb = ylim self.xlog, self.ylog = xlog, ylog def X(self, v): a, b = self.xa, self.xb if self.xlog: v, a, b = math.log10(v), math.log10(a), math.log10(b) return self.x0 + (v - a) / (b - a) * self.w def Y(self, v): a, b = self.ya, self.yb if self.ylog: v = math.log10(max(v, 1e-300)) a, b = math.log10(a), math.log10(b) return self.y0 + self.h - (v - a) / (b - a) * self.h def txt(x, y, s, anchor="middle", size=11, cls="lab"): return ('%s' % (cls, x, y, anchor, size, s)) def poly(ax, pts, cls): s = " ".join("%.1f,%.1f" % (ax.X(a), ax.Y(b)) for a, b in pts) return '' % (cls, s) def frame(ax, xticks, yticks, xfmt="%g", yfmt="%g", xlab="", ylab="", title=""): out = [] for t in xticks: out.append('' % (ax.X(t), ax.y0, ax.X(t), ax.y0 + ax.h)) for t in yticks: out.append('' % (ax.x0, ax.Y(t), ax.x0 + ax.w, ax.Y(t))) out.append('' % (ax.x0, ax.y0, ax.x0, ax.y0 + ax.h, ax.x0 + ax.w, ax.y0 + ax.h)) for t in xticks: out.append(txt(ax.X(t), ax.y0 + ax.h + 16, xfmt % t)) for t in yticks: out.append(txt(ax.x0 - 8, ax.Y(t) + 4, yfmt % t, anchor="end")) if xlab: out.append(txt(ax.x0 + ax.w / 2.0, ax.y0 + ax.h + 34, xlab, size=12)) if ylab: out.append('%s' % (ax.x0 - 44, ax.y0 + ax.h / 2.0, ax.x0 - 44, ax.y0 + ax.h / 2.0, ylab)) if title: out.append(txt(ax.x0, ax.y0 - 12, title, anchor="start", size=12, cls="ttl")) return out def write_figures(outdir, st, analytic_half, d_star, mc, peak_rows): mean, sd, se, qs, conv, zed = mc if not os.path.isdir(outdir): os.makedirs(outdir) # ---- FIGURE 1: the three routes and their residual ----------------------- c, ka = 0.100, 10.0 ** (-PKA_ACETIC) buf = [(c, ka)] ax = Ax(74, 30, 500, 170, (0, 14), (1e-4, 3.0), ylog=True) g = frame(ax, [0, 2, 4, 6, 8, 10, 12, 14], [1e-4, 1e-3, 1e-2, 1e-1, 1e0], xfmt="%g", yfmt="%.0e", xlab="pH", ylab="beta (mol/L per pH)", title="0.100 M acetic acid, pKa 4.756, 25 C") pts = [] p = 0.0 while p <= 14.0001: pts.append((p, beta_multi(10.0 ** (-p), buf, KW_25))) p += 0.05 g.append(poly(ax, pts, "sA")) for (pp, be, bn, bs, rn, rs) in FIG["validation"]: g.append('' % (ax.X(pp), ax.Y(bs))) g.append('' % (ax.X(PKA_ACETIC), ax.y0, ax.X(PKA_ACETIC), ax.y0 + ax.h)) g.append(txt(ax.X(PKA_ACETIC) + 6, ax.y0 + 14, "pKa", anchor="start", size=11, cls="lab2")) g.append(txt(ax.X(1.4), ax.Y(0.6), "water wall", anchor="middle", size=11, cls="lab2")) g.append(txt(ax.X(12.6), ax.Y(0.6), "water wall", anchor="middle", size=11, cls="lab2")) ax2 = Ax(74, 246, 500, 70, (0, 14), (1e-14, 1e-6), ylog=True) g += frame(ax2, [0, 2, 4, 6, 8, 10, 12, 14], [1e-13, 1e-10, 1e-7], xfmt="%g", yfmt="%.0e", xlab="pH", ylab="|rel. diff|") fd = [(pp, max(abs(rn), 1e-14)) for (pp, be, bn, bs, rn, rs) in FIG["validation"]] ts = [(pp, max(abs(rs), 1e-14)) for (pp, be, bn, bs, rn, rs) in FIG["validation"]] g.append(poly(ax2, fd, "sC")) g.append(poly(ax2, ts, "sB")) for pp, v in ts: g.append('' % (ax2.X(pp), ax2.Y(v))) g.append(txt(ax2.X(7.0), ax2.Y(3e-8), "titration route", size=11, cls="lab2")) g.append(txt(ax2.X(7.0), ax2.Y(4e-13), "finite difference", size=11, cls="lab2")) _dump(outdir, "fig1.svg", g, 640, 362) # ---- FIGURE 2: the family, the window ---------------------------------- ax = Ax(74, 30, 500, 250, (0, 14), (0, 0.075)) g = frame(ax, [0, 2, 4, 6, 8, 10, 12, 14], [0.0, 0.015, 0.03, 0.045, 0.06, 0.075], xfmt="%g", yfmt="%.3f", xlab="pH", ylab="beta (mol/L per pH)", title="0.100 M monoprotic buffers, 25 C") # half-capacity band on the pKa 7 buffer b7 = [(0.100, 1.0e-7)] pm = find_beta_max(b7, KW_25, 7.0) bm = beta_multi(10.0 ** (-pm), b7, KW_25) lo, hi = half_capacity_window(b7, KW_25, pm, bm) g.append('' % (ax.X(lo), ax.y0, ax.X(hi) - ax.X(lo), ax.h)) cls = ["sA", "sB", "sC", "sD", "sA", "sB", "sC"] for i, pka in enumerate((3.0, 5.0, 7.0, 9.0, 11.0)): bb = [(0.100, 10.0 ** (-pka))] pts = [] p = 0.0 while p <= 14.0001: v = beta_multi(10.0 ** (-p), bb, KW_25) pts.append((p, min(v, 0.075))) p += 0.06 g.append(poly(ax, pts, cls[i])) g.append(txt(ax.X(pka), ax.Y(0.0605), "%.0f" % pka, size=11, cls="lab2")) g.append('' % (ax.X(lo), ax.Y(bm / 2), ax.X(hi), ax.Y(bm / 2))) g.append(txt(ax.X(7.0), ax.Y(bm / 2) - 8, "half-capacity window %.3f pH wide" % (hi - lo), size=11, cls="lab")) g.append(txt(ax.X(7.0), ax.Y(0.0682), "peak ln10*C/4 = %.5f" % (LN10 * 0.1 / 4), size=11, cls="lab2")) _dump(outdir, "fig2.svg", g, 640, 330) # ---- FIGURE 3: two-component ------------------------------------------ ax = Ax(74, 30, 500, 250, (4.0, 10.0), (0, 0.065)) g = frame(ax, [4, 5, 6, 7, 8, 9, 10], [0.0, 0.013, 0.026, 0.039, 0.052, 0.065], xfmt="%g", yfmt="%.3f", xlab="pH", ylab="beta (mol/L per pH)", title="0.100 M total, split between two buffers about pH 7") for i, (d, cl) in enumerate([(0.0, "sA"), (d_star, "sB"), (2.0, "sC"), (3.0, "sD")]): bb = [(0.050, 10.0 ** (-(7.0 - d / 2))), (0.050, 10.0 ** (-(7.0 + d / 2)))] pts = [] p = 4.0 while p <= 10.0001: pts.append((p, beta_multi(10.0 ** (-p), bb, KW_25))) p += 0.04 g.append(poly(ax, pts, cl)) top = beta_multi(10.0 ** (-7.0), bb, KW_25) g.append(txt(ax.X(7.0), ax.Y(top) - 7, "D = %.3f" % d, size=11, cls="lab2")) g.append('' % (ax.X(7.0), ax.y0, ax.X(7.0), ax.y0 + ax.h)) _dump(outdir, "fig3.svg", g, 640, 330) # ---- FIGURE 4: open vs closed ----------------------------------------- ax = Ax(80, 30, 494, 250, (6.0, 8.4), (1e-3, 1e0), ylog=True) g = frame(ax, [6.0, 6.4, 6.8, 7.2, 7.6, 8.0, 8.4], [1e-3, 1e-2, 1e-1, 1e0], xfmt="%.1f", yfmt="%.0e", xlab="pH", ylab="beta (mol/L per pH)", title="carbon dioxide system, 37 C, pCO2 40 mmHg") op = [(p, bo) for (p, bo, bc, hc) in FIG["carb"]] cl_ = [(p, bc) for (p, bo, bc, hc) in FIG["carb"]] g.append(poly(ax, op, "sA")) g.append(poly(ax, cl_, "sC")) g.append('' % (ax.X(7.4), ax.y0, ax.X(7.4), ax.y0 + ax.h)) g.append('' % (ax.X(7.4), ax.Y(st["open"]))) g.append('' % (ax.X(7.4), ax.Y(st["closed"]))) g.append(txt(ax.X(6.10), ax.Y(0.030), "open, pCO2 held", anchor="start", size=12, cls="lab")) g.append(txt(ax.X(6.10), ax.Y(0.00135), "closed, CT held", anchor="start", size=12, cls="lab")) g.append(txt(ax.X(7.43), ax.Y(0.0012), "blood pH 7.40", anchor="start", size=11, cls="lab2")) g.append(txt(ax.X(7.55), ax.Y(0.0130), "x %.1f" % st["ratio"], anchor="middle", size=12, cls="lab")) g.append('' % (ax.X(7.4), ax.Y(st["open"]), ax.X(7.4), ax.Y(st["closed"]))) _dump(outdir, "fig4.svg", g, 640, 330) # ---- FIGURE 5: Monte Carlo convergence --------------------------------- lo_y = min(m - 2.2 * s_ for (_, m, s_) in conv[3:]) hi_y = max(m + 2.2 * s_ for (_, m, s_) in conv[3:]) ax = Ax(80, 30, 494, 250, (10, 200000), (lo_y, hi_y), xlog=True) g = frame(ax, [10, 100, 1000, 10000, 100000], [round(lo_y + (hi_y - lo_y) * k / 4.0, 3) for k in range(5)], xfmt="%d", yfmt="%.3f", xlab="trials", ylab="open / closed ratio", title="Monte Carlo running mean, seed %d" % MASTER_SEED) up = [(i, m + s_) for (i, m, s_) in conv] dn = [(i, m - s_) for (i, m, s_) in conv] band = " ".join("%.1f,%.1f" % (ax.X(a), ax.Y(b)) for a, b in up) band += " " + " ".join("%.1f,%.1f" % (ax.X(a), ax.Y(b)) for a, b in reversed(dn)) g.append('' % band) g.append(poly(ax, [(i, m) for (i, m, s_) in conv], "sA")) g.append('' % (ax.x0, ax.Y(mean), ax.x0 + ax.w, ax.Y(mean))) g.append(txt(ax.X(300), ax.Y(mean) - 8, "final mean %.4f" % mean, anchor="start", size=11, cls="lab")) g.append(txt(ax.X(20000), ax.Y(lo_y + (hi_y - lo_y) * 0.12), "+/- 1 standard error", anchor="middle", size=11, cls="lab2")) _dump(outdir, "fig5.svg", g, 640, 330) sys.stderr.write("figures written to %s\n" % outdir) def _dump(outdir, name, parts, w, h): with open(os.path.join(outdir, name), "w", encoding="utf-8") as f: f.write('\n') f.write("\n".join(parts)) f.write("\n\n") # ============================================================================== # MAIN # ============================================================================== def main(): svgdir = None if "--svg" in sys.argv: svgdir = sys.argv[sys.argv.index("--svg") + 1] section_header() section_sanity() section_validate_derivative() peak_rows = section_maximum() analytic_half, win_rows, retained = section_window() section_concentration() d_star, two_rows, real_rows = section_two_component() st = section_carbonate() section_no_optimum(st) section_dose(st) mc = section_monte_carlo() section_sensitivity(st) section_summary(st, analytic_half, d_star, mc) if svgdir: write_figures(svgdir, st, analytic_half, d_star, mc, peak_rows) if __name__ == "__main__": main()