#!/usr/bin/env python3 """PD20 conservation-by-construction integrator (the rebuild's step engine). Owns a CLOSED UNIVERSE of pools (serum, peritoneal layers, tissue, apparatus circuit, permeate, plus a `source` supply and a `sink` waste). Every physical process is a paired transfer from physics_core, so the total water and total per-solute mmol across ALL pools are invariant by construction (Cori is the one species conversion, tracked explicitly). step() evaluates every flux from a TOP-OF-STEP SNAPSHOT, then applies them (explicit Euler: gradients at t, state to t+dt), then asserts conservation and halts on any violation. Leaves the legacy simulator.py untouched. Controller/thirst/geometry/pretrain plug into these pools later — they consume the physics, they don't change it. Run: python3 physics_sim.py (closed-system + external-balance + apparatus tests) """ import math, os from physics_geometry import (cavity_axes, height_at_volume, slab_geometry, lumen_pressure_cmH2O, volume_from_pressure_L) from physics_core import (SOLUTES, REFLECTION, SIEVING, R_T, Pool, VALENCE, transfer_fluid, transfer_solute, total_water, total_solute, ro_nominal_perm, conductivity, cp_beta, osmotic_pressure_diff_conc, MACRO_OSM_FACTOR, PROTEINS, PROTEIN_MW_G_PER_MMOL, RO_PRESSURE_MIN_BAR, RO_PRESSURE_MAX_BAR) # ---- whole-body solute model constants (added 2026-08-02, advisory review; see working_on.md) K_IC_PER_KG = 0.04 / 70.0 # intercompartmental K clearance, L/min PER KG. Weight-scaled 2026-08-02 # after an adversarial pass: a fixed 0.04 L/min against a weight-scaled serum # volume made the redistribution constant 3.0 h at 45 kg and 7.0 h at 90 kg — # backwards. Same for the capacity term below. K_IC = 0.04 # (legacy scalar, kept for any caller that reads it; step() uses K_IC_PER_KG) # HONEST PROVENANCE (corrected 2026-08-02 after adversarial review): this is NOT # identified by the Fernandez acute-load test as originally claimed here. That # test's peak spans the whole 1.06 +- 0.13 band for K_IC anywhere in [0.02, 0.15], # and is INSENSITIVE to S (peak 5.147/5.148/5.147 at S = 100/250/600). What # actually reproduces the ESRD peak is VD_FRACTION['Potassium'] = 0.20. K_IC and S # are therefore ORDER-OF-MAGNITUDE choices inside a physiologic range, not fitted # constants: S from the clinical 200-400 mEq-per-mEq/L rule, K_IC set so # redistribution runs in hours. Treat both as sensitivity parameters. K_BODY_PER_PLASMA_PER_KG = 250.0 / 70.0 # S per kg. Weight-scaled 2026-08-02: unscaled, the 45 kg # anuric who DEFINES the potassium constraint was over-buffered 14.3x versus the # 90 kg patient's 7.4x. K_BODY_PER_PLASMA = 250.0 # S: mmol of body K per mmol/L of plasma K. CAPACITY of the buffer. SOURCE: INFERRED/UNVERIFIED (order-of-magnitude buffer capacity) # The clinical rule "200-400 mEq total body deficit per 1 mEq/L fall in serum K" # (NCBI Clinical Methods, 'Serum Potassium'). Documented as non-linear there # (~200-300 at 1 mEq/L, ~500-600 at 2), so a constant S is a first-order fit. # ---- SKELETAL (bone) reservoir for the mineral species, added 2026-08-04 ----------------------- # Bone holds ~99% of body calcium, ~85% of phosphate and ~60% of magnesium, and it is a genuine # exchangeable buffer on a days-to-weeks timescale (rapid surface exchange, not just remodelling). # Without it, mineral intake has nowhere to go: adding sourced dietary Ca/Mg drove serum ionised # calcium to 2.36 (ref 1.15-1.30) because dialysate removal was the ONLY exit. Same structural trap # the potassium intracellular pool was hiding until 2026-08-02. # Form is the SAME linearised buffer used for potassium: the equilibrium set point rides on bone # CONTENT, so at rest J == 0 exactly (a buffer, never a source or a chronic drain), and as the pool # depletes the set point falls and the buffer weakens. # C_eq = C0 + (Q_bone - Q_bone0)/S ; J = K_BONE * (C_eq - C_ecf) [+ve = bone -> serum] # Body contents (mmol, 70 kg): Ca ~25000, PO4 ~19000, Mg ~1000. S values are the body-store change # per unit serum change; these are ORDER-OF-MAGNITUDE choices inside physiologic ranges, NOT fitted # constants — the mineral literature gives no clean intercompartmental clearance, exactly as it # gives none for potassium. Treat as sensitivity parameters and label them as such. BONE_CONTENT = {'Calcium': 25000.0, 'Phosphate': 19000.0, 'Magnesium': 1000.0} # mmol at 70 kg BONE_K = {'Calcium': 0.60, 'Phosphate': 0.40, 'Magnesium': 0.30} # L/min per 70 kg BONE_S = {'Calcium': 400.0, 'Phosphate': 300.0, 'Magnesium': 150.0} # mmol per mmol/L # ---- NON-OSMOTIC SODIUM STORE (skin + skeletal-muscle GAG-bound Na), added 2026-08-24 (S19/P2) -- # 23Na-MRI and balance studies (Titze J; Kopp C; Dahlmann A) show a substantial fraction of body Na # sits in skin and muscle interstitium bound to negatively-charged glycosaminoglycans, osmotically # INACTIVE: it accumulates without commensurate water. The model carried body Na = ECF only # (tissue['Sodium'] = 0 — measured 1960 mmol), so 100% of every retained mmol expanded ECF # osmolality, which is the structural half of the D6 hypernatraemia artifact. # Same linearised buffer as bone/K, zero water carried (osmotically silent by construction): # C_eq = Na0 + (Q_skin - Q_skin0) / (S_Na * wt/70); J = K_Na * (C_eq - C_serumNa) [+ve = skin->serum] # K_Na = serum.volume / K_DAYS (a time constant of K_DAYS days; audit: 3-7). # Content = S_Na * Na0 (~1000 mmol at 70 kg — the audit's "missing 1000-1500 mmol of buffering"), # which buffers a fraction S_Na/(S_Na + Vd_Na) = 7/21 ~ 33% of an acute Na load osmotically # inactive (audit's 20-40%), Vd_Na = 14 L. Both knobs overridable so every sodium claim is # reportable across the sweep; PD20_SKIN_NA_S=0 disables the store. SKIN_NA_S = float(os.environ.get('PD20_SKIN_NA_S', '7.0')) # mmol per mmol/L, scaled by wt/70 SKIN_NA_K_DAYS = float(os.environ.get('PD20_SKIN_NA_K_DAYS', '5.0')) # ECF exchange time constant, days ALB_SYN_HALFLIFE_MIN = 1.0 * 1440.0 # restorative albumin synthesis time constant (1 d) SOURCE: INFERRED/UNVERIFIED (net reserve above baseline turnover) ALB_SYN_MAX_MMOL_MIN = 10.0 / 66500.0 * 1000.0 / 1440.0 # 10 g/d NET-synthesis reserve, MW 66.5 kDa. SOURCE: INFERRED/UNVERIFIED (10 g/d net reserve) # The reserve is the EXTRA synthesis above baseline turnover available to # compensate peritoneal protein loss. Measured effluent albumin loss is ~6.5 g/d, # so a 3 g/d reserve left serum albumin declining ~20% over 30 days; 10 g/d holds # it near baseline (the physiologic hepatic maximum is ~15 g/d). VD_FRACTION = {'Sodium': 0.20, 'Chloride': 0.20, 'Glucose': 0.20, 'Potassium': 0.20, 'Urea': 0.60, 'Albumin': 0.20, 'Globulin': 0.20, 'Icodextrin': 0.20, 'Bicarbonate': 0.42, 'Lactate': 0.20, 'Calcium': 0.20, 'Magnesium': 0.20, 'Phosphate': 0.30, 'Sulfate': 0.25, 'Creatinine': 0.60} # total-body-water distributed, like urea (lit ~0.6 L/kg) # Phosphate 0.30: ECF plus the rapidly-exchangeable fraction. LIMITATION, stated: the # large SKELETAL and intracellular phosphate reservoirs are NOT modelled, so serum # phosphate here is MORE responsive to intake and removal than in a real patient. Do # not read absolute phosphate control performance off this without that caveat — it is # the same trap the potassium Vd was hiding until 2026-08-02. # Ca/Mg: only the IONISED extracellular fraction is modelled, so ECF (0.20). The vast # skeletal reservoir is NOT modelled — these are short-horizon balance species here. # HONEST PROVENANCE (annotated 2026-08-22, literature audit finding #4). Chagnac's measurement # constrains anatomic CONTACT AREA against fill volume; this exponent is then applied to MTAC, which # is a different quantity. Chagnac's own creatinine mass-transfer figures quoted nearby (112 -> 142 # mg for +46% volume) would imply ~0.627 — essentially the geometric 2/3 this value was chosen # over — but a larger dialysate volume removes more mass at CONSTANT MTAC, so that figure is not an # MTAC scaling either. NO PUBLISHED SOURCE PINS THE MTAC EXPONENT. The two candidates differ by # 1.46^0.19 ~ 7.6% in MTAC at a 3 L fill, so this is a real but bounded uncertainty. # STATUS: INFERRED / UNVERIFIED — do not present it in any publication as measured for MTAC. AREA_RECRUIT_EXP = 0.437 # exponent in effective exchange area proportional to fill_volume^n. # MEASURED, not geometric. Chagnac et al. measured peritoneal surface area # CONTACT (PSA-CD) directly at two fill volumes in PD patients: going from a # 2 L to a 3 L instil raised intraperitoneal volume 46 +/- 2% and increased # PSA-CD by only 18 +/- 2.3% (0.57 -> 0.67 m2), with creatinine mass transfer # rising 112 -> 142 mg. Solving 1.46^n = 1.18 gives n = 0.437. # The engine previously used the GEOMETRIC exponent 2/3, which predicts +29% # area for that same +46% volume — a 1.6x overestimate of how much membrane a # larger fill actually recruits. The geometric value assumes the cavity is a # smooth expanding vessel; real peritoneum is already largely in contact at # clinical fills, and the extra fluid mostly deepens the existing pool rather # than wetting proportionally more membrane. # THIS MATTERS BEYOND ACCURACY: the fill-volume lever is claimed as an # area-recruitment mechanism for small-solute clearance, and 2/3 overstated # its benefit. Any conclusion resting on raising fill volume must be re-read. # INTRAPERITONEAL PRESSURE vs FILL, measured in adult CAPD patients lying flat. # Durand 1993 (PMID 8105960): 13 cmH2O at 2.82 L. Durand 1994 (PMID 7999866): linear in volume from # 2 to 5 L, 18 cmH2O the maximal acceptable pressure. Two points on that line — 13 at 2.82 L and the # ~4 cmH2O/L slope the two papers imply together — give the intercept below. IAP_SLOPE = 4.0 # cmH2O per litre of intraperitoneal volume (Durand 1994, linear to 5 L) IAP_BASE = 13.0 - IAP_SLOPE * 2.82 # = 1.72 cmH2O at zero fill (Durand 1993 anchor) PERIT_MIN_WET_L = 0.15 # cavity volume (L) below which peritoneal transport is scaled toward zero. SOURCE: INFERRED/UNVERIFIED (residual film; ~250 mL per Oberg & Rippe) # A drained cavity has no dialysate to exchange with, and real peritoneum # retains a residual film after any drain. See the wetted-volume guard in # step(). Inert above this volume, i.e. in all continuous-flow operation. PERIT_MIX_D = 1.2e-6 # effective inter-layer diffusivity in the peritoneal cavity (m^2/s). SOURCE: INFERRED/UNVERIFIED (calibrated to a 10-min mixing time) # Lumps CONVECTIVE mixing (respiration, posture, peristalsis) with molecular # diffusion (glucose ~6.7e-10 m^2/s, so this is ~1.8e3 x molecular). # CALIBRATION (corrected 2026-08-07): the slowest mode of a diffusing column # of height h decays with tau = h^2/(pi^2 * D). A 70 kg cavity at a 2.45 L # fill stands h = 0.085 m, so a 10-minute cavity mixing time requires # D = h^2/(pi^2 * 600 s) = 1.2e-6 m^2/s. The previous 1.5e-5 was 12x too # large and gave a mixing time under one minute — it was never detected # because an explicit-flux overshoot cap, not this constant, was setting the # actual rate. SENSITIVITY PARAMETER, not a measured constant — no study # reports an in-vivo peritoneal eddy diffusivity. Raising it collapses the # cavity to one well-mixed pool; lowering it stratifies. Swept in # test_layer_mixing.py. DWELL_SLEW = 0.15 # max peritoneal fill/drain rate (L/min = 150 mL/min) — physiologic pump speed; # slews the tidal so dwell changes are smooth (no step at meal onset / window exit). # ACTUATOR ENVELOPES (inventor 2026-08-15: "no changes ever happen instantaneously, they take time, # so speed increased and decreased linearly, same for glucose concentration, mixer addition, # pressure, all has envelopes"). Until now DWELL_SLEW above was the ONLY rate-limited actuator: the # dialysate FLOW stepped instantly between baseline and post-meal boost, and the mixer's dextrose and # sodium, the RO pressure and the recovery setpoint were all applied in full on the step they were # commanded. An actuator with no envelope can traverse its whole range in ONE step, so shrinking the # control tick makes the modelled device strictly FASTER than any real one — precisely the wrong # direction for a claim, and it gets worse as the tick gets finer. # # THESE RATES ARE ENGINEERING ASSUMPTIONS, NOT MEASUREMENTS. No pump or mixer datasheet has been # consulted. Each is chosen so the actuator crosses its full working range in a few minutes, which is # the right order for the hardware class; every one is overridable so the inventor can drop in real # numbers without touching this file. Set PD20_SLEW=off to restore the old instantaneous behaviour # (which is what every number published before today was measured with). _SLEW_ON = os.environ.get('PD20_SLEW', '').strip().lower() not in ('off', '0', 'false', 'no') def _slew_rate(name, default): raw = os.environ.get('PD20_SLEW_' + name.upper(), '') if not raw.strip(): return float(default) v = float(raw) if v <= 0.0: raise ValueError('PD20_SLEW_%s must be > 0 (use PD20_SLEW=off to disable envelopes ' 'entirely), got %r' % (name.upper(), raw)) return v # RATES ONLY — NO ENDPOINTS, AND NO TRAVERSAL TIMES. # These comments twice leaked controller set-points out of a file that is published verbatim. # 1st pass (2026-08-20): the dextrose CAP and the sodium BAND endpoints were named outright. # Removed — but the removal was incomplete, in two ways an adversarial review then caught: # 2nd pass (2026-08-20, same day): # (a) the FLOW line three rows up still named "0.005-0.070", which IS FLOW_MIN/FLOW_MAX from # physics_nn_controller.py — a WITHHELD module's constants, published verbatim. # (b) removing the numerals but keeping a traversal time left the ARITHMETIC intact: rate x # time recovers the endpoint to within a couple of percent. A redaction that leaves the # product is not a redaction. The same applied to the RO pump line, which named its target. # (c) 2026-08-22: the comment written to explain (b) ITSELF SPELLED OUT the multiplication and # named the constant, so the publish gate refused physics_sim.py — correctly. A note about # a leak is still the leak if it contains the numbers. Describe the CLASS, never the values. # So: state the RATE, which is physics (how fast an actuator can move), and nothing about the LIMIT # it travels toward, which is a controller set-point. No endpoints, no spans, no traversal times. FLOW_SLEW = _slew_rate('flow', 0.020) # dialysate flow, L/min per min DEX_SLEW = _slew_rate('dex', 20.0) # mixer dextrose, mmol/L per min SOURCE: design choice (actuator slew) NA_SLEW = _slew_rate('na', 5.0) # mixer sodium, mmol/L per min PRESSURE_SLEW = _slew_rate('pressure', 20.0) # RO pump, bar per min RECOVERY_SLEW = _slew_rate('recovery', 0.05) # reject-valve position, fraction per min def _ramp(target, current, rate, dt): """Move `current` toward `target` at no more than `rate` per minute over `dt` minutes. Linear both directions, which is what the inventor specified ("speed increased and decreased linearly"). current=None means the actuator starts already at its commanded value, so a run does not spend its first minutes ramping up from zero and contaminating the comparison.""" if current is None: return target step = rate * dt return min(current + step, max(current - step, target)) # Per-solute renal concentrating factor applied to the filtered-fraction (gfr_frac·serum) urine # concentration. A residual nephron CONCENTRATES urea far above plasma (residual urea clearance ~ GFR), # so urine urea >> serum; Na/Cl/HCO3/Lactate stay near the filtered fraction (≈1); glucose ~0 # (reabsorbed). Potassium is cleared by the distal-secretion term (renal_k), so it keeps ×1 here. URINE_CONC_MULT = {'Urea': 10.0, 'Glucose': 0.05, 'Albumin': 0.0, 'Globulin': 0.0, # SOURCE: INFERRED/UNVERIFIED (per-solute urine concentration multiplier) # Divalents and phosphate are NOT handled like Na/Cl. Calcium is ~98-99% reabsorbed (fractional # excretion 1-2%), magnesium ~95% (FE 3-5%), while PHOSPHATE is the opposite — the kidney is the # main phosphate exit and fractional excretion RISES steeply as GFR falls, exceeding 1.0 relative to # the filtered fraction in advanced CKD. Flat multipliers of 1.0 were wrong for all three. 'Calcium': 0.10, 'Magnesium': 0.30, 'Phosphate': 3.0, 'Sulfate': 3.0} # Albumin 0.0: a residual nephron at ~3 mL/min filters essentially no albumin (66.5 kDa). The # generic urine route was carrying it at serum concentration, losing 13.2 g/d — more than twice the # physiologic PERITONEAL loss. Stool albumin is likewise zero (see stool_conc below). Together # these two spurious routes were 19.7 of the 27.2 g/d total albumin loss. Added 2026-08-02. # Cell-membrane reflection σ for the serum<->tissue (ECF<->ICF) water shift. Unlike # the leaky peritoneal membrane (σ_Na≈0.05), the cell membrane is near-impermeant to # Na/K/glucose (σ≈1, they drive tonicity), while urea crosses freely (σ≈0, no # tonicity effect). H3 fix — was incorrectly reusing the peritoneal σ=0.05. CELL_SIGMA = {'Sodium': 1.0, 'Chloride': 1.0, 'Glucose': 1.0, 'Potassium': 1.0, # SOURCE: dead code (superseded) 'Urea': 0.0, 'Bicarbonate': 1.0, 'Lactate': 0.5} # (superseded by the algebraic # transcellular water PARTITION; kept for reference) # Insulin-driven transcellular potassium shift (Na-K-ATPase), as a self-returning BUFFER that # tracks its OWN displacement (an amount, `ik_displaced`), NOT ICF concentration — so it is # fully decoupled from the water partition (which changes tissue volume, hence [K]icf). Net # serum->ICF = gain*insulin_pert*[K]serum − K_RETURN*ik_displaced. First term: a glucose RISE # above fasting pumps K into cells, scaled by insulin_sensitivity (diabetic shifts LESS). # Second term: relaxes the displacement back to zero (stable under sustained hyperglycemia; # self-zeroing when glucose returns to fasting). Zero at fasting for ANY patient -> no chronic # dump, no 1/si blowup, no water-partition coupling. K_INS_GAIN = 0.02 # insulin-K pump gain (per min) SOURCE: INFERRED/UNVERIFIED (order-of-magnitude) K_RETURN = 0.002 # displacement return rate (per min); sets the transient magnitude SOURCE: INFERRED/UNVERIFIED (order-of-magnitude) class _TrueVolumeView: """Concentration view for a pool whose solutes really are dissolved in its own water (the ICF compartment). conc = solutes / pool.volume, no apparent-Vd scaling.""" def __init__(self, pool): self.pool = pool def __getitem__(self, s): return self.pool.volume class _VdView: """Per-solute distribution volume that SCALES with the compartment's current water volume (ECF=0.20·BW reference). So Vd_s = pool.volume·(vd_frac_s/0.20); concentration = solutes/Vd then responds correctly to water gain/loss (matching the legacy per-solute Vd model, but kept consistent with the live volume instead of frozen). Conservation is unaffected (it's mmol-based).""" def __init__(self, pool): self.pool = pool def __getitem__(self, s): return self.pool.volume * (VD_FRACTION[s] / 0.20) # --------------------------------------------------------------------------------------------- # CIRCUIT SELECTION BY ENVIRONMENT (added 2026-08-13). # The closed-loop architecture has to be trainable and testable end to end, and the sim is # constructed in dozens of places — the rollout generator, the cohort runner, every experiment. # Threading a configuration argument through all of them would touch far more code than the change # is worth and would silently miss any callsite that was overlooked, which is exactly how a # controller ends up TRAINED on one circuit and MEASURED on another. One environment switch, # applied at the end of construction, reaches every callsite by definition. # UNSET = nothing happens and the engine is bit-for-bit what it was. # PD20_CIRCUIT=closed full recirculation + 2 L reservoir + multi-species bed # PD20_CIRCUIT=ro the shipped reverse-osmosis circuit (explicit, same as unset) # Individual overrides, applied after the preset: PD20_REGEN, PD20_RESERVOIR_L, PD20_BLEED_LD, # PD20_SORBENT_MULTI, PD20_RO_PASSES. def env_ro_passes(): """PD20_RO_PASSES as an int >= 1, or None when unset or blank. VALIDATED, NOT COERCED. `int(os.environ.get(...) or default)` accepted "0" — a truthy string — and produced ro_passes=0, which then behaved EXACTLY like 1 because the recycler loop runs `range(max(0, ro_passes - 1))`. A configuration silently becoming a different configuration is the same defect class as the two this module already had (a pressure that did not recalibrate the membrane, and a pass count that did not reach the run), so this one fails loudly instead. Blank is treated as unset because that is what an exported-but-empty shell variable means.""" raw = os.environ.get('PD20_RO_PASSES') if raw is None or not raw.strip(): return None try: n = int(raw.strip()) except ValueError: raise ValueError('PD20_RO_PASSES must be an integer >= 1, got %r' % raw) if n < 1: raise ValueError('PD20_RO_PASSES must be an integer >= 1 (the RO stage always runs at ' 'least one pass), got %r' % raw) return n VALID_REGEN = ('ro', 'sorbent', 'ro+sorbent', 'sorbent+ro') def _validate_regen(value): """`regen` must be one of the four modes the step engine actually branches on.""" v = str(value).strip() if v not in VALID_REGEN: raise ValueError('regen must be one of %s, got %r. A near-miss silently degrades to plain ' 'reverse osmosis with the sorbent stage absent, which is a different ' 'device.' % (', '.join(VALID_REGEN), value)) return v def _apply_env_circuit(sim): # Applied pressure is collected here and committed ONCE at the end through set_ro_pressure(), # never assigned directly. ro_sp (per-solute permeance) and ro_ref_jw (reference water flux) are # both back-solved FROM the pressure, so a bare assignment leaves the membrane calibrated for one # pressure while being driven at another: at 80 bar on a 15-bar calibration the realised urea # rejection is 0.911 against the nominal 0.500, which scrubs the recirculated dialysate and # flatters clearance, and ff = jw_raw/ro_ref_jw pins at 1.0 so it never throttles. None = the # pressure was never touched, so nothing is recomputed and the default engine is untouched. want_pressure = None preset = os.environ.get('PD20_CIRCUIT', '').strip().lower() if preset == 'closed': sim.regen = 'sorbent' sim.reservoir_cap_L = 2.0 sim.sorbent_multi = True elif preset == 'highrec': # THE CONFIGURATION THAT WORKS. Reverse osmosis KEPT — it is the loop's only exit for the # osmotic agent, and removing it turns the circuit into a one-way glucose trap — but run at # high recovery so little water leaves with the concentrate, with the serial sorbent column # closing the urea gap the membrane leaves and the reservoir recovering the tidal drain. sim.regen = 'ro+sorbent' sim.reservoir_cap_L = 2.0 want_pressure = 80.0 elif preset in ('ro', 'open'): sim.regen = 'ro' sim.reservoir_cap_L = 0.0 sim.sorbent_multi = False _g = os.environ.get if _g('PD20_REGEN'): sim.regen = _validate_regen(_g('PD20_REGEN')) if _g('PD20_RESERVOIR_L'): sim.reservoir_cap_L = float(_g('PD20_RESERVOIR_L')) if _g('PD20_BLEED_LD'): sim.loop_bleed_L_day = float(_g('PD20_BLEED_LD')) if _g('PD20_SORBENT_MULTI'): sim.sorbent_multi = _g('PD20_SORBENT_MULTI') not in ('0', 'false', 'False', '') if _g('PD20_TIDAL_TO_LOOP'): sim.tidal_to_loop = _g('PD20_TIDAL_TO_LOOP') not in ('0', 'false', 'False') _passes = env_ro_passes() if _passes is not None: sim.ro_passes = _passes if _g('PD20_RO_PRESSURE'): want_pressure = float(_g('PD20_RO_PRESSURE')) if want_pressure is not None: sim.set_ro_pressure(want_pressure) # recalibrates ro_sp AND ro_ref_jw (PF6) class PhysicsSim: def __init__(self, weight=70.0, N=10, enable=None, mtac_scale=1.0, sex='M', fill_L=None, circuit_holdup_L=0.0, regen='ro'): self.N = N self.weight = weight self.mtac_scale = mtac_scale # transporter type: >1 high, <1 low self.sex = sex self.circuit_holdup = circuit_holdup_L # device priming/dead volume (L); 0 = instantaneous self.enable = {'peritoneal': True, 'apparatus': True, 'external': True, 'cori': True} if enable: self.enable.update(enable) # SEX-based total body water (Watson): male TBW ~0.60·wt, female ~0.50·wt. ECF:ICF held # at 1:2 (ECF = 1/3 TBW, ICF = 2/3 TBW). Male reproduces the prior 0.20/0.40 split EXACTLY # (regression-preserving default). (FULL_PHYSICS_AUDIT §G.) tbw_frac = 0.60 if str(sex).upper().startswith('M') else 0.50 ecf = (1.0 / 3.0) * tbw_frac * weight # serum/ECF water (L) tissue0 = (2.0 / 3.0) * tbw_frac * weight # ICF water (L) # PERITONEAL FILL VOLUME is now a parameter so HIGH/LOW fill can be simulated # (FULL_PHYSICS_AUDIT §F). Default 0.035·wt (~2.45 L @ 70 kg). Recruited membrane area # (hence MTAC and UF coefficient Kf) scales with fill ∝ volume^AREA_RECRUIT_EXP (measured, # referenced to the default fill so default behavior is UNCHANGED (factor = 1 at default). ref_fill = 0.035 * weight perit0 = fill_L if fill_L is not None else ref_fill # peritoneal fill (L) self.fill_L = perit0 self.fill_area_factor = (perit0 / ref_fill) ** AREA_RECRUIT_EXP if ref_fill > 0 else 1.0 # Physiologic starting concentrations (mmol/L). Not tuned to any outcome. serum_c = {'Sodium': 140, 'Chloride': 103, 'Glucose': 5, 'Potassium': 4, 'Urea': 20, 'Albumin': 0.7, 'Globulin': 0.187, 'Icodextrin': 0, 'Bicarbonate': 24, 'Lactate': 1, 'Calcium': 1.20, 'Magnesium': 0.85, 'Phosphate': 1.60, 'Sulfate': 2.1, # lumped uraemic anion (sulfate + organic acids), divalent. 'Creatinine': 0.70} # ~8 mg/dL, ESRD-typical (normal 0.6-1.3 mg/dL = 53-115 umol/L, Wikipedia) # ELECTRONEUTRAL AT t=0 (fixed 2026-08-22, S19/P5). Sulfate was 5.0 mmol/L = 10 mEq/L, # standing in for the anion gap WHILE albumin also carried its full Figge charge # (0.7 mmol/L x 18.6 = 13.0 mEq/L) — the two double-counted, so the initial serum was # -5.8 mEq/L (anion excess) and drifted +2.5 mEq/L/day. The charge balance # Na+K+2Ca+2Mg = Cl+HCO3+Lac+1.8*PO4+2*SO4+18.6*Alb fixes sulfate at # (148.1 - 143.9)/2 = 2.1 mmol/L. The lumped solute still represents the real # unmeasured-anion pool (sulfate ~0.3-0.5 mmol/L plus organic acids ~3-4 mEq/L), # just at the charge-correct total instead of a double-counted one. # ESRD serum phosphate 1.6 mmol/L (~5.0 mg/dL) — typical pre-dialysis hyperphosphataemia # serum IONISED Ca 1.15-1.30 and Mg 0.7-1.0 mmol/L (clinical reference ranges) # Initial peritoneal fill — ELECTRONEUTRAL: Cl = Na+K−HCO3−Lactate = 132+2−35−10 = 89 # (was 97, an 8 mEq anion excess). Flushed out quickly by the apparatus, but kept balanced. dial_c = {'Sodium': 132, 'Chloride': 89, 'Glucose': 0, 'Potassium': 2, 'Urea': 0, 'Albumin': 0, 'Globulin': 0, 'Icodextrin': 0, 'Bicarbonate': 35, 'Lactate': 10, 'Calcium': 1.25, 'Magnesium': 0.25, 'Phosphate': 0.0, 'Sulfate': 0.0, 'Creatinine': 0.0} # dialysate carries none — it is a solute to REMOVE # dialysate phosphate is ZERO — phosphate is a solute to REMOVE, never supplied # Physioneal-class divalents # init total-pool mmol = conc · Vd, with Vd = volume·(vd_frac/0.20) self.serum = Pool(ecf, {s: serum_c[s] * ecf * (VD_FRACTION[s] / 0.20) for s in SOLUTES}) # ICF solute content. Only POTASSIUM is a real intracellular reservoir in this engine (set # below). Every other solute was initialised at serum_c * V * (VD_FRACTION/0.20), which put # 19.6 mmol of ALBUMIN inside cells (albumin is strictly extracellular) and 3x the # physiologic exchangeable sodium (5880 vs ~2000-2500 mmol) into a pool nothing reads. # Zeroed 2026-08-02: the tissue concentrations were computed into an unused local and # discarded, and the osmotic partition uses the separate `icf_osmoles` scalar, so this is # behaviourally inert — verified against test_published_regression. self.tissue = Pool(tissue0, {s: 0.0 for s in SOLUTES}) # Per-solute distribution volumes that scale with current pool volume. self.serum_vd = _VdView(self.serum) # The ICF pool IS the compartment its solutes occupy, so its concentration view must be its # own water volume — NOT _VdView, which scales by VD_FRACTION/0.20 and is only meaningful # for SERUM (an ECF pool carrying an apparent whole-body Vd). Applied to the ICF pool it # reported potassium at 70 mmol/L in a 28 L compartment holding 3920 mmol, i.e. exactly half # the physiologic 140. Found by advisory review 2026-08-02; pinned by test_icf_pool_sane.py. self.tissue_vd = _TrueVolumeView(self.tissue) # ICF (tissue) as a REAL intracellular compartment: high-K reservoir (~140 mM direct, # ~3500 mmol) for the insulin-K pump, plus a tracked effective-osmole content for the # transcellular water PARTITION. Init ISO-TONIC to serum so the resting state has zero # water shift (the partition places serum/ICF at osmotic equilibrium each step). self.tissue.solutes['Potassium'] = 140.0 * tissue0 # ICF K reservoir self.icf_osmoles = (2.0 * serum_c['Sodium'] + serum_c['Glucose']) * tissue0 # iso-tonic baseline self.ik_displaced = 0.0 # insulin-K cumulative displacement (mmol) self.layers = [Pool(perit0 / N, {s: dial_c[s] * (perit0 / N) for s in SOLUTES}) for _ in range(N)] # Circuit is the device dead/priming volume (well-mixed holdup). Primed with dialysate # so it is at steady standing volume from t=0; holdup=0 -> Pool(0.0) = prior behavior. self.circuit = Pool(circuit_holdup_L, {s: dial_c[s] * circuit_holdup_L for s in SOLUTES}) self.permeate = Pool(0.0) # Large external reservoirs so the universe is closed and finite. self.source = Pool(1e6, {s: 1e6 * serum_c[s] for s in SOLUTES}) # intake+enrichment supply self.sink = Pool(0.0) # urine+insensible+reject # HOLDING RESERVOIR (added 2026-08-13). Capacity 0 = ABSENT, which is exactly the prior # circuit, so every published number is unchanged unless this is switched on deliberately. # It exists because waste attribution showed the dominant loss is not the membrane: the # tidal dwell lever drains the belly straight to the drain when the setpoint falls, and # refills it from fresh supply when the setpoint rises, several times a day. The fluid is # thrown away and immediately re-bought. A reservoir closes that cycle — it stores the # regenerated surplus while the belly is small and gives it back when the belly grows. self.reservoir = Pool(0.0) self.reservoir_cap_L = 0.0 # TWO SEPARATE DECISIONS, PREVIOUSLY WELDED TOGETHER. Recovering the tidal shrink-drain is # (a) a ROUTING choice — send the surplus into the circuit to be regenerated instead of down # the drain, which is a valve, not a vessel — and (b) a STORAGE choice — hold regenerated # fluid the belly cannot accept yet, which does need a tank. The engine gated (a) on # `reservoir_cap_L > 0`, so the only way to get the routing was to accept a 2 L tank. # None = legacy behaviour (routing follows the tank); True/False sets the valve explicitly, # so a tankless circuit that still recovers its own tidal drain can be measured. self.tidal_to_loop = None # PURGE / BLEED LINE (added 2026-08-13, L/day; 0 = absent = prior circuit exactly). # Measured reason it exists: with the loop fully closed, DEXTROSE HAS NO EXIT. The mixer can # only ADD solute, so once the recirculating fluid carries more glucose than the target # recipe the controller has no authority to bring it down — and it cannot even see the # problem, being concentration-blind. probe_reservoir_death.py caught this cleanly: the # commanded dextrose was IDENTICAL to three decimals with and without the reservoir # (14.514 vs 14.517) while dialysate glucose ran from 81 to 93 mmol/L, ultrafiltration ran # away and the patient died of fluid depletion on day 6.9. # Under reverse osmosis the reject stream is that exit: rejecting glucose at 0.99 is not # only a water loss, it is the loop's only glucose bleed. Closing the loop removes it. # A purge valve restores it explicitly and at a chosen, auditable cost in litres. self.loop_bleed_L_day = 0.0 self.sorbent_multi = False # recirculating bed removes urea only (True = full per-species) self.V_dwell = perit0 # physiological peritoneal dwell volume the self.V_dwell_avg = perit0 # running average of V_dwell (for control_tbw compensation # under continuous tidal — the instantaneous V_dwell swings # with the tidal cycle, but the controller should see the # TIME-AVERAGED device-instilled excess, not the per-step # value, or it over-ultrafiltrates when the belly drains low. # apparatus holds (real PD keeps the belly filled; # net fluid removal = osmotic UF drained out, NOT # cavity drainage). Replaces the old empty-belly # hydrostatic-UF artifact + the artificial net_removal knob. # Net device-commanded residual fill(+)/drain(-) rate (L/min) = J_in - J_out, i.e. the # controller-commanded TIDE (LUK-004 claims 41-49): >0 => outflow rate < inflow rate # (post-prandial net fill); <0 => outflow > inflow (return to baseline). Set each step # from the slew-limited dwell change; its time-integral equals the residual-volume rise. self.resid_dVdt_Lmin = 0.0 # RO membrane calibration (Wijmans-Baker), Kf, lymph, MTAC, Fickian self.ro_wp = 5.0 self.ro_pressure = 15.0 # Concentration-polarisation modulus at the reference flux. 1.0 = OFF, and OFF is an exact # no-op that reproduces every published number byte-for-byte (asserted by # test_polarization.py). Real RO runs ~1.1-1.4; at >=1.1 the 2-pass configuration stops # meeting the concentration-blind mixer bound (exp_polarization.py), so this is opt-in and # changing it is an inventor decision, not a default. self.ro_cp_beta = 1.0 self.ro_sp, self.ro_ref_jw = ro_nominal_perm(self.ro_wp, ro_pressure=self.ro_pressure) # RECYCLER passes: how many times the recovered fluid goes through the RO membrane # (LUK-003 claim 27 recirculation loop). 1 = single pass (prior behavior). >1 re-filters # the permeate to strip more waste solute (urea) from the re-infused dialysate -> better # clearance, at the cost of recovering less water per unit effluent. self.ro_passes = 1 self.ro_recovery_capped = False # set each step: did the osmotic wall bind? # Dialysate regeneration mode: 'ro' (RO split: recover clean water, DUMP the # concentrate incl. glucose) or 'sorbent' (REDY/wearable-kidney sorbent: strip # urea to waste, RETAIN glucose, recirculate — glucose-sparing regeneration). # VALIDATED, NOT COERCED (2026-08-18, adversarial review). `regen` was compared by string # equality in three places with no membership check and no else-branch, so 'ro_sorbent', # 'RO+SORBENT' or 'banana' all silently became plain RO with the sorbent stage ABSENT — # measured: sorbent urea bound 6.63 mmol on 'ro+sorbent' vs 0.000 on each typo. Two of the # valid values ('ro+sorbent' and 'sorbent+ro') are near-identical strings with OPPOSITE # placement semantics and a 4.6x difference in cartridge loading, so a typo is both easy and # consequential. Same defect class env_ro_passes() was written to eliminate, on the adjacent # knob. Fail loudly instead. self.regen = _validate_regen(regen) self.sorbent_urea_strip = 0.95 # sorbent urea-removal fraction per pass # SERIAL SORBENT POLISHING COLUMN (regen='ro+sorbent', added 2026-08-12). # PLACEMENT AND WHY: the column sits DOWNSTREAM of the reverse-osmosis stage and treats the # PERMEATE, not the raw effluent. That is the efficient placement for one specific reason — # reverse osmosis is very good at the salts and very POOR at urea. This engine's own nominal # rejections are Na/Cl ~0.995 (commercial polyamide, see physics_core.ro_nominal_perm), # macromolecules >=0.99, but urea only 0.50, and that 0.50 # is real membrane physics rather than a modelling shortcut: a small neutral polar molecule # of 60 Da passes a polyamide membrane readily. So half the urea the device recovers is # re-infused into the patient, and urea is the very solute dialysis exists to remove. # A sorbent column placed on the permeate therefore attacks exactly the gap the RO leaves, # and it does so on the SMALL clean stream rather than the whole effluent, which is what # keeps the cartridge small enough to be body-worn. # PER-SPECIES REMOVAL, from the classical multi-layer sorbent train (urease converting urea # to ammonium carbonate, a cation exchanger binding the ammonium, a hydrous-oxide layer for # phosphate and other anions, and activated carbon for the middle-molecule organics): self.sorbent_strip = {'Urea': 0.95, 'Sulfate': 0.60, 'Phosphate': 0.55, 'Potassium': 0.30} # CARTRIDGE LOADING COUNTER (mmol bound, per species). A sorbent bed is consumed by the MASS # it binds, not by how long it runs or how fast fluid goes past, so this is the quantity that # decides cartridge life and it is the only way to answer whether slowing the flow preserves # the cartridge. Accumulated wherever the column strips, in either placement. self.cum_sorbent_bound = {s_: 0.0 for s_ in SOLUTES} # Species the column must NOT take, because they are prescribed and would have to be # re-dosed downstream: glucose, sodium, chloride, bicarbonate, lactate, calcium, magnesium, # albumin and globulin all pass through untouched. GLUCOSE RETENTION is the point of a # sorbent stage — a carbon bed does not adsorb glucose appreciably, so unlike the reject # stream the permeate keeps its osmotic agent. # LIMITATION, STATED: a real cation-exchange layer RELEASES sodium (and hydrogen) in # exchange for the ammonium and potassium it binds, which is the best-known drawback of # sorbent regeneration. That exchange is NOT modelled here, so this arm OVERSTATES the # benefit to the extent that a real column would add a sodium load the mixer must then # subtract. Treat the numbers as an upper bound until the exchange stoichiometry is added. # TOTAL PERITONEAL FLUID ABSORPTION — renamed 2026-08-09; it is NOT lymph flow. # The value 1.07 mL/min is correct and measured, but the NAME was wrong by a factor of 3-5: # Rippe B, Rosengren BI, Venturoli D, Microcirculation 2001;8(5):303-20 (PMID 11687943) puts # TRUE lymphatic absorption of macromolecules at "approximately 0.2 mL/min", and Oberg & # Rippe 2017 (PMC5733752, Table 1) list "Peritoneal lymph flow (L) (ml/min): 0.3", against a # TOTAL fluid absorption clearance of ~1 mL/min. What this term models — bulk fluid leaving # the cavity by every route — is the ~1 mL/min total, of which lymph proper is a minority # and direct tissue uptake the rest. The attribute keeps the name `lymph_rate` for now # because renaming it touches every caller; the PHYSICAL QUANTITY is total absorption and # any publication text must say so. # Was 0.10 L/h = 1.67 mL/min, which is not a literature value and sits ABOVE the measured # resting range. Imholz AL, Koomen GC, Struijk DG, Arisz L, Krediet RT, "Effect of an # increased intraperitoneal pressure on fluid and solute transport during CAPD", Kidney Int # 1993;44(5):1078-85 (PMID 8264138), measured the lymphatic # absorption rate DIRECTLY in CAPD patients: 1.07 +/- 0.18 mL/min at rest, rising to # 1.86 +/- 0.25 mL/min under abdominal compression. The resting figure is the quantity this # constant represents, so it is used. # WHY IT MATTERED: lymphatic absorption is a near-constant drain on the cavity, so an # overestimate subtracts the same volume from every exchange regardless of bag strength. # A single-exchange validation against published PET kinetics showed the engine # under-producing ultrafiltration by a roughly constant 300-450 mL per 4 h dwell at every # dextrose strength — the exact signature of a rate error rather than a broken mechanism. # KNOWN REMAINING GAP: the same study shows this rate is PRESSURE-DEPENDENT (+74% under # compression). This model holds it constant, so it does not capture the penalty that a # larger intraperitoneal fill imposes on net ultrafiltration. That coupling is relevant to # any fill-volume conclusion and is not yet implemented. self.lymph_rate = 1.07 / 1000.0 # L/min (1.07 mL/min, PMID 8264138) # DYNAMIC membrane area: MTAC (small solutes) and Kf recruit with the CURRENT dwell VOLUME # (area ∝ V^AREA_RECRUIT_EXP), recomputed each step from V_dwell — raising the peritoneal fill # increases K/urea/UF clearance in real time (the clearance lever flow can't be, since K is # MTAC/area-limited, not flow-limited). Base (per-unit-area) coefficients kept so the effective # values track V_dwell. Macromolecules (albumin/icodextrin) do not scale. At a CONSTANT dwell = # initial fill this reproduces the prior frozen values EXACTLY (regression-safe). self.ref_fill = ref_fill # SODIUM and CHLORIDE 0.010 -> 0.0045 L/min on 2026-08-09. Oberg CM & Rippe B, Kidney Int # Rep 2017;2(6):1128-38 (PMC5733752, Table 1) give a DIRECT human value: "PS for Na+ and # anions (ml/min): 4.5". The engine ran 10.0 mL/min, 2.2x that, justified only by an # INFERENCE ("bracketed by creatinine 9.6 and urea 22.9") which was never a measurement of # sodium at all. POTASSIUM is left at 0.010 and is now labelled INFERRED, not verified: no # reachable source gives a peritoneal PS or MTAC for potassium, so pretending otherwise # would be the same mistake in the opposite direction. # GLUCOSE MTAC 0.012 -> 0.0090 L/min (9.0 mL/min), calibrated 2026-08-22 against the # Twardowski PET by test_pet_twardowski.py. This is RECORD.md open item O11 and it is now # closed on the reference bag. # WHY IT MOVED: at 12 mL/min the membrane absorbed glucose far too fast — D/D0 at 4 h read # 0.31 against a published 0.38 for 2.5%, the osmotic agent was gone before the dwell ended, # and net UF came out 225 mL short. 9.0 mL/min sits inside the published human range # (~8-14 mL/min by MTAC methods; the old 12 was near the top of it) and brings the 2.5% # reference bag into tolerance on BOTH axes at once: D/D0 0.391 vs 0.38, net UF +280 vs +400 # (tol 150). 4.25% UF also lands, +842 vs +800. # WHY THIS IS CALIBRATION AND NOT FITTING: exactly ONE parameter moved, and it moved the two # INDEPENDENT published observables — solute equilibration and water movement — toward their # targets TOGETHER. A fitted parameter buys one at the other's expense; measured across the # sweep, D/D0 and UF improved monotonically in step. RECORD.md's standard is met: "changing # two parameters to land inside a range is fitting, not calibrating." # NOT AN OPEN DEFECT — an earlier version of this comment claimed one, wrongly, and the # correction is kept here because the error is easy to make twice. It read: "D/D0 is nearly # flat across bag strength (0.41/0.39/0.37) where published falls steeply (0.49/0.38/0.26)". # That comparison is INVALID. 0.49/0.38/0.26 are Twardowski's TRANSPORTER CATEGORIES — low, # average, high — all measured with the SAME 2.5% bag, which is the only strength the # standard PET uses. They are not D/D0 by bag strength. D/D0 really is only weakly dependent # on strength, so flat is CORRECT; the strength axis is where UF varies, and it does. # Measured against the right scale, the model matches the category curve within tolerance: # 0.500 / 0.391 / 0.294 vs 0.49 / 0.38 / 0.26. See test_pet_twardowski.py, which tests the # two axes separately for exactly this reason. self._mtac_unit = {'Sodium': 0.0045, 'Chloride': 0.0045, 'Glucose': 0.0090, 'Potassium': 0.010, 'Urea': 0.020, 'Bicarbonate': 0.010, 'Lactate': 0.010, 'Calcium': 0.006, 'Magnesium': 0.006, 'Creatinine': 0.010, # 10 mL/min, the D/P-creatinine PET standard (Krediet/Amsterdam) # Phosphate clears POORLY in PD (large hydrated radius, protein binding): # real CAPD removes only ~300-400 mg/d against a ~700 mg/d absorbed load, # which is why hyperphosphataemia persists and binders are needed. # PHOSPHATE 0.005 -> 0.0102 L/min (fixed 2026-08-22, literature audit # finding #3). 0.005 was a factor of TWO below the PS 10.2 mL/min in the # SAME published table this module quotes for sodium, and it was justified # from a clinical-OUTCOME argument ("real CAPD removes only ~300-400 mg/d") # rather than from a transport measurement. Fitting a transport coefficient # to an outcome hides whatever else is producing that outcome — here the # 0.30 distribution volume and the un-modelled skeletal/intracellular # phosphate pool, both already self-flagged. Use the measured PS; let the # outcome fall out of the physics. 'Phosphate': 0.0102, 'Sulfate': 0.0102} # protein-bound uraemic anions clear poorly self._Kf_unit = 8e-5 self.mtac = {'Albumin': 0.0001, 'Globulin': 0.00005, 'Icodextrin': 0.0001} self._apply_dwell_area() # sets self.mtac (small solutes) from V_dwell # Kf FIXED at the initial-fill area (NOT re-scaled by later dwell changes — UF stays # controller-managed and undisturbed by the dwell clearance lever). self.Kf = self._Kf_unit * self.fill_area_factor # NOTE: effective membrane surface is encoded in MTAC (glucose 9.0 mL/min, urea 20 mL/min — # lit ~8-14/15-25) and Kf (0.08 mL/min/mmHg — lit 0.03-0.09), validated by # test_pet_twardowski.py against the published PET on both axes. (Glucose read "12 mL/min" # here for a while after the constant above was calibrated to 9.0 — a stale duplicate of a # number that lives three lines up. Restating a constant in prose is how that happens.) The old explicit geometric area_cm2=1000 (0.1 m², # too small vs real ~1-2 m²) drove only the inter-layer Fickian, which the well-mixed cavity # replaced — removed as dead code. # Conservation bookkeeping over the closed universe. self.sim_min = 0.0 # wall-clock (min since start) for meal timing / dwell schedule self.meal_events = [] # absolute sim_min of ad-hoc "I ate" events (GUI); [] = clock-only self.cori_consumed = 0.0 # Lactate removed by Cori self.cori_produced = 0.0 # HCO3 added by Cori # Clinical tracking (pure bookkeeping, conservation-neutral): # Conductivity readings — the ONLY device sensing of the fluid (bulk, not concentrations). # Measured each step at three points: effluent drained from the belly, recovered permeate # before the mixer, and the regenerated dialysate infused back. The controller sees these # (and their lag statistics), never any serum concentration. self.cond_out = 0.0 # drained peritoneal effluent conductivity (mS/cm) self.cond_pre = 0.0 # recovered permeate, pre-mixer (post-RO, pre-enrichment) self.cond_in = 0.0 # regenerated dialysate infused into the belly self.cum_uf = 0.0 # net water removed from serum across the peritoneal # membrane (osmotic UF − back-filtration − lymph), L. # +ve = fluid ultrafiltered off the patient. self.cum_glu_abs = 0.0 # net glucose absorbed serum<-dialysate across the # membrane (diffusion + sieved convection), mmol. # WHERE THE DISCARDED WATER ACTUALLY GOES (added 2026-08-13). The sink volume was a single # number, so "the device wastes N L/day" could not be attributed to a mechanism — and I # mis-read that number three times in one session. There are exactly three volume paths into # the sink and they have completely different engineering answers, so they are counted apart: # RO REJECT the concentrate the membrane could not recover. Irreducible for a given # recovery, and ZERO under full recirculation (regen='sorbent'). # EXCESS DUMP fully regenerated dialysate discarded because the cavity had no room for # it at that instant — the dwell setpoint fell, or ultrafiltration had # already refilled the belly. This is NOT a membrane loss: the fluid was # clean and ready to infuse. A holding reservoir could return it later. # REMAINDER numerical leftovers, must stay at floating-point noise. # Sorbent stripping moves SOLUTE only and adds no volume, so it appears in none of these. # ACTUATOR STATE for the envelopes. None = "not yet commanded", so the FIRST step adopts the # commanded value outright and the ramp applies only to CHANGES thereafter. self._act_flow = None self._act_dex = None self._act_na = None self._act_recovery = None self._act_enrich = {} # the dialysate recipe ACTUALLY infused this step self.cum_reject_L = 0.0 # RO concentrate discarded, L # UREA REMOVAL ATTRIBUTED PER STAGE (mmol, cumulative). "How much clearance does the FIRST # RO pass do, how much does the SECOND, and how much does the sorbent column do?" cannot be # answered from totals — each stage is in series on the same stream, so the only honest # split is to book the urea each one actually sends to waste at the moment it does it. # The sorbent's share is already in cum_sorbent_bound['Urea']. These add the rest. self.cum_urea_ro1 = 0.0 # urea leaving in the FIRST RO pass's reject stream self.cum_urea_ro2p = 0.0 # urea leaving in the reject of passes 2..n (the recycler) self.cum_urea_bleed = 0.0 # urea leaving in the deliberate loop purge self.cum_urea_excess = 0.0 # urea leaving with regenerated fluid the cavity refused self.cum_urea_remainder = 0.0 # urea in the numerical remainder swept to waste # TIDAL SHRINK-DRAIN was an UNATTRIBUTED waste path: the per-stage counters above summed to # 8.4% less urea than actually reached the sink over 24 h of tidal cycling, and the gap was # booked nowhere. Any "which stage does the clearance" answer was short by that much. self.cum_urea_dwell_drain = 0.0 # urea leaving with cavity fluid the falling dwell discards self.cum_excess_dump_L = 0.0 # regenerated dialysate the cavity could not accept, L self.cum_remainder_L = 0.0 # numerical remainder swept to waste, L # TIDAL SHRINK-DRAIN. When the controller LOWERS the dwell setpoint, the surplus cavity # fluid is removed — and it goes STRAIGHT TO THE DRAIN, bypassing the regenerator entirely. # Found 2026-08-13 by attribution: it is 14.0 L/day and it is IDENTICAL in every # architecture, because it happens upstream of the regeneration choice. Under reverse # osmosis it hides inside a 65 L/day bill and looks like nothing; under full recirculation # it is 77% of the entire waste stream. Set `reservoir_cap_L` > 0 to route it into the loop. self.cum_dwell_drain_L = 0.0 # cavity fluid discarded by a FALLING dwell setpoint, L self.cum_patient_loss_L = 0.0 # urine + insensible + stool water leaving the patient, L self.cum_reservoir_stored_L = 0.0 # regenerated dialysate parked in the reservoir, L self.cum_reservoir_returned_L = 0.0 # reservoir fluid given back to the patient, L self.cum_loop_bleed_L = 0.0 # deliberate purge of the recirculating loop, L self.cum_drain_deficit_L = 0.0 # commanded drain the cavity could NOT supply (L, cumulative) self.drain_deficit_steps = 0 # steps on which that happened — 0 in all normal operation self.cum_settle_residual_L = 0.0 # layer-settling residual left when the pass budget ran out self.settle_unconverged_steps = 0 # NOTE: these four counters are diagnostics. They are asserted zero in normal operation by # test_layer_mixing.py and are NOT currently surfaced in run_cohort / report output — an # unsatisfiable drain is recorded but would still be invisible in a production run. self.cum_glu_drained = 0.0 # glucose DRAINED from the peritoneal cavity through the outflow # pump (to circuit→RO→waste), mmol — glucose that LEFT the belly # but was NOT absorbed by the patient. The PATIENT only absorbs # what crosses the membrane; what drains out is wasted. For the # glucose economics: added = absorbed + drained_waste + perit_remaining. self.cum_glu_rejected = 0.0 # glucose LEAVING the loop to WASTE (RO reject + excess- # permeate dump), mmol — the drained dialysate glucose NOT # recovered. Bookkeeping only (conservation-neutral). self.cum_glu_added = 0.0 # FRESH glucose drawn from supply (source) into the # dialysate (enrich + makeup), mmol. THE patent metric: # recycling makes this << single-pass Dianeal. self.cum_na_added = 0.0 # FRESH sodium drawn from supply into the dialysate # (enrich + makeup), mmol — dialysate THROUGHPUT, NOT what # the patient receives (most leaves with the effluent). self.cum_cl_added = 0.0 # FRESH chloride the MIXER draws from supply (enrich + makeup), # mmol — mixer consumption tracking (per-day/30-day grams reported). self.cum_drink_L = 0.0 # oral water intake, L (thirst-driven; a RESULT, not an input) self.cum_na_diet = 0.0 # dietary sodium into the patient (oral intake), mmol. # NON-OSMOTIC SODIUM STORE (S19/P2) — a real Pool with zero water, registered in _pools() # so conservation is structural, exactly like bone. Content scales with weight; the resting # content is S_Na * Na0 (~1000 mmol at 70 kg). self._na0 = self.serum.conc('Sodium', self.serum_vd['Sodium']) # resting serum Na (C0) self.skin = Pool(0.0, {'Sodium': SKIN_NA_S * self._na0 * weight / 70.0}) self._skin0 = self.skin.solutes['Sodium'] # Q_skin0 (reference) self._na_rate_hist = [] # (sim_min, serum Na) samples for the 24h Na-rate fatal guard self._patient_na0 = (self.serum.solutes['Sodium'] + self.tissue.solutes['Sodium'] + self.skin.solutes['Sodium']) # body Na @ t0 (incl. store) # skeletal pool, scaled to body weight (bone mass scales with size) # Skeletal reservoir as a REAL Pool with zero water volume, registered in _pools(). # It was a plain dict until 2026-08-05; that made conservation depend on every caller # REMEMBERING to add it, and test_icf.py (which legitimately rebaselines _S0 after # injecting solute) did not — producing a spurious 'Calcium non-conservation +2.5e4 mmol'. # As a Pool it is seen by total_solute()/total_water() automatically and every bone<->serum # movement goes through the paired transfer_solute primitive, so conservation is structural # rather than remembered. # OVOID cavity geometry (added 2026-08-07). The cavity is an ellipsoid filled from the # floor, so fluid HEIGHT is non-linear in volume and horizontal slabs of equal height hold # very unequal volumes — 51 mL at the floor vs 497 mL at the top of a 3.4 L fill, with the # area-to-volume ratio spanning 11x. That is what lets a real inter-layer gradient exist. self.cav = cavity_axes(weight) self.overfill_L = 0.0 # peritoneal volume beyond the ovoid's capacity (L) self._geo_cache = (None, None, None) # (rounded volume, h, slabs) — see _geometry() # PD_Z_INFLOW / self.z_inflow REMOVED 2026-08-18 (adversarial review). Three defects on one # line and zero consumers: repo-wide grep found exactly ONE reference — the assignment # itself — so the comment "None -> top of fluid" described behaviour nothing implemented; # `or None` swallowed the legitimate value 0.0; and an exported-but-EMPTY PD_Z_INFLOW raised # ValueError and killed every PhysicsSim() construction, the exact case env_ro_passes() # documents and handles correctly a few lines away. A dead switch that looks live is worse # than no switch: it invites an experiment that silently does nothing. self.z_outflow = 0.0 # drain lumen sits on the cavity floor (gravity drainage) self.bone = Pool(0.0, {sp: BONE_CONTENT[sp] * weight / 70.0 for sp in BONE_CONTENT}) self._bone0 = {sp: self.bone.solutes[sp] for sp in BONE_CONTENT} self._min0 = {sp: serum_c[sp] for sp in BONE_CONTENT} self._alb0 = self.serum.solutes['Albumin'] # restorative albumin-synthesis set point self._glb0 = self.serum.solutes['Globulin'] # and the same for globulin (see 5h-bis) self._k0_plasma = serum_c['Potassium'] # resting plasma K (C0) for the ICF K buffer self._q_icf0 = self.tissue.solutes['Potassium'] # resting ICF K content (Q_icf0) self._W0 = total_water(self._pools()) self._S0 = {s: total_solute(self._pools(), s) for s in SOLUTES} # PER-SOLUTE t0 BASELINE for net_retained(). Placed HERE, not beside _patient_na0: # bone is constructed later than that line, so capturing there raised AttributeError. self._patient0 = {_s: (self.serum.solutes.get(_s, 0.0) + self.tissue.solutes.get(_s, 0.0) + self.bone.solutes.get(_s, 0.0) + self.skin.solutes.get(_s, 0.0)) for _s in SOLUTES} self.tbw0 = self.control_tbw() # baseline as the CONTROLLER estimates it (kept for the # controller's own target and for callers that report the # device's view of fluid balance) self.patient_tbw0 = self.patient_tbw() # baseline TRUE patient water (cavity excluded) — # the basis for the volume-overload fatal check. Must not be # read through control_tbw(): its V_dwell_avg EMA cannot # track a cycler, so the trip fired on estimator lag rather # than on the patient (adversary round 6, S2b). _apply_env_circuit(self) def _apply_dwell_area(self): """Recompute MTAC (small solutes) from the CURRENT dwell volume: recruited diffusive peritoneal surface area ∝ V^AREA_RECRUIT_EXP (Chagnac, measured), referenced to the default fill. Kf (the UF coefficient) recruits by the SAME law — LpS and PS are the same area parameter in the three-pore model, so a dwell change moves both. (A constant dwell = init is unchanged.)""" self.fill_area_factor = (self.V_dwell / self.ref_fill) ** AREA_RECRUIT_EXP if self.ref_fill > 0 else 1.0 a = self.mtac_scale * self.fill_area_factor for s in ('Sodium', 'Chloride', 'Glucose', 'Potassium', 'Urea', 'Bicarbonate', 'Lactate', 'Calcium', 'Magnesium', 'Phosphate', 'Sulfate', 'Creatinine'): self.mtac[s] = self._mtac_unit[s] * a # Kf RECRUITS WITH THE SAME AREA (fixed 2026-08-22, literature audit finding #2). # Kf was deliberately frozen at the initial fill while MTAC recruited, described as a # "clean separation: dwell = clearance lever, dextrose = UF lever". That separation is not # physical. In the three-pore model LpS and PS are THE SAME area parameter (A0/dx), and the # Chagnac contact-area argument this very method invokes for MTAC applies identically to # LpS. As written, raising the dwell volume bought clearance at ZERO ultrafiltration cost — # which no membrane does — so every conclusion about the fill-volume lever inherited a free # lunch. Both now scale by fill_area_factor, from one law. self.Kf = self._Kf_unit * self.mtac_scale * self.fill_area_factor def net_na_retained(self): """Net sodium the PATIENT actually gained (body = serum + ICF + non-osmotic store), mmol since t0. ~0 at steady state = balanced (dietary intake matched by device removal); this is the real 'salt given to the patient', NOT the dialysate throughput (cum_na_added). The skin store is included because it is PATIENT sodium, just osmotically inactive.""" return (self.serum.solutes['Sodium'] + self.tissue.solutes['Sodium'] + self.skin.solutes['Sodium'] - self._patient_na0) def net_retained(self, species): """Net amount of `species` the PATIENT actually gained (body = serum + ICF + bone), mmol since t0. The generalisation of net_na_retained to any solute. Added 2026-08-22 on the inventor's standing reporting rule: every run must report total gained salt, gained Na, gained Cl and gained glucose PER DAY, alongside clearance. Only sodium had a net-retention accessor; chloride and the rest were reported only as MIXER THROUGHPUT (`cum_cl_added`), which is what the DEVICE circulated, not what the PATIENT gained — two quantities that differ by orders of magnitude and must never be confused in a publication. Bone is included because Ca/Mg/PO4 partition into it. """ b = self.bone.solutes.get(species, 0.0) sk = self.skin.solutes.get(species, 0.0) return (self.serum.solutes.get(species, 0.0) + self.tissue.solutes.get(species, 0.0) + b + sk - self._patient0.get(species, 0.0)) def _pools(self): return [self.serum, self.tissue, *self.layers, self.circuit, self.permeate, self.source, self.sink, self.bone, self.skin, self.reservoir] def ro_max_recovery(self, feed, beta=1.0): """Highest fraction of the feed that reverse osmosis can recover before the concentrate's own osmotic pressure equals the applied pressure — the OSMOTIC WALL. Extracting water leaves the solute behind, so at recovery r the concentrate is 1/(1-r) times the feed concentration and its osmotic pressure is pi_feed/(1-r). Water stops crossing when that reaches the applied pressure: pi_feed / (1 - r_max) = P_applied => r_max = 1 - pi_feed / P_applied WHY THIS WAS NEEDED (2026-08-13). `perm_vol = recovery * ff * feed.volume` took the COMMANDED recovery — 0.7 — with no reference to the concentrate at all. Measured on the live engine: spent dialysate 9.12 bar against 15 bar applied, so the true ceiling is 39.2%, and at 70% the concentrate would have needed 30.4 bar against 15 applied — a net driving pressure of MINUS 15.4 bar. The engine was recovering water against a gradient that should have stopped it. The pre-existing `ff` (flux-feasible) factor does NOT catch this: it evaluates the flux at the FEED concentration against a reference flux, so it never sees the concentrate. The wall is a property of how far you concentrate, not of where you start — which is exactly why a cap written on the feed missed it. beta is the concentration-polarisation modulus: the membrane WALL sees more than the bulk, so polarisation pulls the ceiling in further.""" osm = 0.0258 * sum(feed.conc(_s) * MACRO_OSM_FACTOR.get(_s, 1.0) for _s in SOLUTES) if self.ro_pressure <= 0.0: return 0.0 return max(0.0, 1.0 - beta * osm / self.ro_pressure) def _perit_volume(self): return sum(L.volume for L in self.layers) def peritoneal_conc(self, solute): """Cavity-mean concentration of a solute (mmol/L). For tracking peak/mean peritoneal glucose — the invention holds this LOW & STEADY where a Dianeal bolus spikes HIGH then decays at equal time-average.""" v = self._perit_volume() return (sum(L.solutes[solute] for L in self.layers) / v) if v > 1e-12 else 0.0 def perit_glu_mmol(self): """Total glucose currently sitting in the peritoneal cavity (mmol). This is glucose that has been instilled but NOT YET absorbed and NOT YET drained — it is still in the belly. THE GLUCOSE ECONOMY (corrected 2026-08-07 — the identity stated here before was wrong): cum_glu_added == cum_glu_abs + cum_glu_rejected + perit_glu_mmol() + loop holdup where loop holdup is the glucose standing in `circuit` + `permeate`. Verified to close at 0.000 mmol over a 48 h run. `cum_glu_drained` is NOT a term in that identity. Drained glucose returns to the cavity via the circuit, so it is a THROUGHPUT counter (how much glucose has cycled out of the belly), not a SINK. The only true sinks are absorption into the patient and RO reject to waste. Including it would double-count. Measured independently: `cum_glu_drained` equals the glucose actually crossing cavity -> circuit to 2.4e-15 relative. The controller is BLIND to all of this (no peritoneal sensor); it is analysis only.""" return sum(L.solutes['Glucose'] for L in self.layers) # ---- the integrator --------------------------------------------------- def set_ro_pressure(self, ro_pressure): """Retune the applied RO pressure AND recompute the flux-feasible reference so the ff = jw_raw/ro_ref_jw cap stays valid (PF6 fix — never set self.ro_pressure directly).""" self.ro_pressure = ro_pressure self.ro_sp, self.ro_ref_jw = ro_nominal_perm(self.ro_wp, ro_pressure=ro_pressure) def step(self, dt, action=None): action = action or {} self.sim_min += dt # advance the wall clock (meal/dwell scheduling) # sodium-rate history for the osmotic-demyelination fatal guard (added 2026-08-24): a 24h # rolling window of serum Na, so fatal_check can bound |dNa|/24h without a second history # mechanism. Trimmed to ~24h so it stays O(1) memory per step. self._na_rate_hist.append((self.sim_min, self.serum.conc('Sodium', self.serum_vd['Sodium']))) while self._na_rate_hist and self._na_rate_hist[0][0] < self.sim_min - 1440.0: self._na_rate_hist.pop(0) # Reset here, not inside the RO branch. The declaration says "set each step", but it was only # assigned where RO actually runs, so a sorbent-only or zero-flow step reported the LAST RO # step's verdict indefinitely. (adversarial review 2026-08-18) self.ro_recovery_capped = False q = action.get('outflow_Lmin', 0.0) # peritoneal drain rate recovery = action.get('recovery', 0.0) # RO recovery setpoint enrich = action.get('enrich', {}) # dialysate electrolyte add (mmol/L of permeate) # RO PRESSURE AS A COMMANDED ACTUATOR (2026-08-22). Absent/None leaves the configured # pressure untouched, so every existing caller, every regression pin and every published # number is byte-identical when nothing commands it. # Applied HERE for the same reason the other envelopes are: step() is the one point the app, # run_nn, the PI teacher and the rollout generator all pass through, so no callsite can be # overlooked and no caller can command past the hardware rating. # Routed through set_ro_pressure() and NEVER by assigning self.ro_pressure directly: the # permeability pair (ro_sp, ro_ref_jw) is recalibrated FROM the pressure, and a bare # assignment silently invalidates the flux-fraction cap (the PF6 defect). _p_cmd = action.get('ro_pressure_bar') if _p_cmd is not None: _p = min(RO_PRESSURE_MAX_BAR, max(RO_PRESSURE_MIN_BAR, float(_p_cmd))) if abs(_p - self.ro_pressure) > 1e-9: self.set_ro_pressure(_p) # RECORD THE RECIPE ACTUALLY INFUSED. Without this there is no way for a test to inspect # what the engine really used — test_actuator_envelopes' electroneutrality check recomputed # chloride from sodium locally and asserted an identity that is algebraically ZERO for any # value, so it passed even with the engine's re-derivation deleted entirely. A guard must # read the engine's output, not its own arithmetic. Set at the END of the envelope block # below, in both modes. # ACTUATOR ENVELOPES — applied HERE, at the one point every caller passes through (the app, # run_nn, the PI teacher and the rollout generator all reach the plant only via step), so no # callsite can be overlooked. Commanded values become SETPOINTS; what the plant actually does # ramps toward them at a bounded rate. See FLOW_SLEW & co. above. # RO PRESSURE as a commandable, envelope-limited actuator. Inert unless a caller passes it # (nothing does today — pressure is set per run), so this changes no existing result; it is # here because the filing discloses the controller adjusting "the pressure used for reverse # osmosis", and a high-pressure pump is the actuator LEAST able to move instantly. _p_cmd = action.get('ro_pressure') if _p_cmd is not None: _p = _ramp(float(_p_cmd), self.ro_pressure, PRESSURE_SLEW, dt) if _SLEW_ON else float(_p_cmd) if abs(_p - self.ro_pressure) > 1e-12: self.set_ro_pressure(_p) # recalibrates ro_sp AND ro_ref_jw (PF6) # The applied-actuator state is recorded EITHER WAY — with envelopes off it simply equals the # command. That keeps `_act_*` meaning "what the plant is actually doing" in both modes, so # instrumentation and tests can read it without first asking which mode is active. if _SLEW_ON: q = self._act_flow = _ramp(q, self._act_flow, FLOW_SLEW, dt) recovery = self._act_recovery = _ramp(recovery, self._act_recovery, RECOVERY_SLEW, dt) else: self._act_flow = q self._act_recovery = recovery if enrich: if _SLEW_ON: _dex = self._act_dex = _ramp(enrich.get('Glucose', 0.0), self._act_dex, DEX_SLEW, dt) _na = self._act_na = _ramp(enrich.get('Sodium', 0.0), self._act_na, NA_SLEW, dt) # ELECTRONEUTRALITY IS NOT OPTIONAL. Ramping sodium without re-deriving chloride would # infuse a dialysate whose cations and anions no longer balance — the one invariant # this engine asserts everywhere. Chloride is recomputed from the SLEWED sodium using # the same identity the mixer used to build the dict (divalents count twice). _na_cmd = enrich.get('Sodium', _na) enrich = dict(enrich) enrich['Glucose'] = _dex enrich['Sodium'] = _na # RE-DERIVE CHLORIDE ONLY IF SODIUM ACTUALLY MOVED (fixed 2026-08-20, adversarial # review). This block was written as the electroneutrality repair for a SLEWED # sodium, but it ran UNCONDITIONALLY whenever envelopes were on — so it silently # replaced the CALLER'S commanded chloride even when nothing was ramping, and only # in that mode. Measured on the engine's own reference recipe with sodium already at # target: commanded Cl 97 became 87 with envelopes on and stayed 97 with them off — # a 10 mmol/L difference in infused fluid, 10.9% in mixer chloride, and 4.8% in # cond_in, which is a LIVE CONTROLLER INPUT. Toggling the envelopes must not change # the dialysate. # When sodium DID move, the re-derivation is required and correct (divalents count # twice). When it did not, the caller's recipe is used as given — and CHECKED, so an # electroneutrality error in a caller is reported instead of silently corrected. if abs(_na - _na_cmd) > 1e-12: enrich['Chloride'] = (_na + enrich.get('Potassium', 0.0) + 2.0 * enrich.get('Calcium', 0.0) + 2.0 * enrich.get('Magnesium', 0.0) - enrich.get('Bicarbonate', 0.0) - enrich.get('Lactate', 0.0)) else: _bal = (_na + enrich.get('Potassium', 0.0) + 2.0 * enrich.get('Calcium', 0.0) + 2.0 * enrich.get('Magnesium', 0.0) - enrich.get('Bicarbonate', 0.0) - enrich.get('Lactate', 0.0)) if abs(_bal - enrich.get('Chloride', _bal)) > 1e-6: raise ValueError( 'dialysate recipe is not electroneutral: Chloride %.4f commanded, %.4f ' 'required by Na+K+2Ca+2Mg-HCO3-Lactate. The engine will not silently ' 'correct a caller recipe.' % (enrich.get('Chloride', float('nan')), _bal)) else: self._act_dex = enrich.get('Glucose', 0.0) self._act_na = enrich.get('Sodium', 0.0) self._act_enrich = dict(enrich) if enrich else {} intake = action.get('intake_Lmin', 0.0) intake_c = action.get('intake_conc', {}) urine = action.get('urine_Lmin', 0.0) gfr_frac = action.get('gfr_frac', 0.6) insensible = action.get('insensible_Lmin', 0.0) cori = action.get('cori_mmol', 0.0) # DWELL-VOLUME lever: the controller may set a new peritoneal fill volume (L). The volume # change is realized PATIENT-NEUTRALLY, moving only DEVICE dialysate: GROW the belly with fresh # dialysate from source, SHRINK it by draining excess to waste. (Doing this via the normal # permeate-recycle refill entangled recovered PATIENT water and dumped it on the shrink, # causing a ~0.75 L over-UF — the bug this fixes.) MTAC then tracks the new volume. self.resid_dVdt_Lmin = 0.0 # reset the commanded tide differential each step dwell_L = action.get('dwell_L', None) if dwell_L is not None: dwell_L = max(0.05, dwell_L) # peritoneal fill is physically > 0 (guard vs NaN area) # SLEW-RATE LIMIT: a real fill/drain pump moves fluid at a bounded rate — cap the dwell # change to <= DWELL_SLEW L/min so the belly RAMPS toward any target instead of stepping. # This makes the tidal genuinely smooth (no base->peak jump at meal onset, no step-drain at # window exit) regardless of what the scheduler requests; ~100 mL/min is a physiologic # instill/drain rate. (Constant dwell = target already met -> no effect.) slew = DWELL_SLEW * dt dwell_L = min(self.V_dwell + slew, max(self.V_dwell - slew, dwell_L)) if dwell_L is not None and abs(dwell_L - self.V_dwell) > 1e-12: dV = dwell_L - self.V_dwell self.resid_dVdt_Lmin = dV / dt # J_in - J_out this step: >0 => outflow < inflow (tide) pv = self._perit_volume() # WHICH LUMEN DOES THE FILL/DRAIN USE? (added 2026-08-08) # PD20 has a double-lumen catheter: it instils through the UPPER lumen and drains through # the LOWER one, which is what creates the vertical gradient. A conventional cycler has a # single access and must both fill and drain through the SAME lumen, and that lumen has # to sit at the bottom of the cavity or the patient cannot be drained. Simulating # conventional APD therefore requires the fill to enter at the FLOOR, not the top — # otherwise the comparator silently inherits PD20's geometry and is not a fair control. # `fill_layer`: None (default) = top lumen, PD20 behaviour, regression-safe. # 0 = bottom lumen, i.e. a conventional single-access catheter. _fl = action.get('fill_layer', None) _tgt = self.layers[self.N - 1] if _fl is None else self.layers[int(_fl)] if dV > 0: # grow belly with fresh dialysate (device fluid) # RESERVOIR FIRST. Stored fluid is already-regenerated dialysate, so the mixer only # has to top up the DEFICIT to the target recipe — the same recycling credit the # permeate path already gets. Fresh supply covers whatever the reservoir cannot. _from_res = min(dV, self.reservoir.volume) if self.reservoir_cap_L > 0 else 0.0 if _from_res > 1e-12: _stage = Pool(0.0) transfer_fluid(self.reservoir, _stage, _from_res) for s in SOLUTES: _def = (enrich.get(s, 0.0) - _stage.conc(s)) * _stage.volume _mvd = min(max(0.0, _def), self.source.solutes[s]) transfer_solute(self.source, _stage, s, _mvd) if s == 'Glucose': self.cum_glu_added += _mvd elif s == 'Sodium': self.cum_na_added += _mvd elif s == 'Chloride': self.cum_cl_added += _mvd transfer_fluid(_stage, _tgt, _stage.volume) self.cum_reservoir_returned_L += _from_res _fresh = dV - _from_res if _fresh > 1e-12: transfer_fluid(self.source, _tgt, _fresh, carried_conc={s: enrich.get(s, 0.0) for s in SOLUTES}) # the tidal FILL draws fresh dextrose/Na/Cl from supply — track it as mixer # consumption (else the FRESH-glucose metric misses the tidal turnover). self.cum_glu_added += _fresh * enrich.get('Glucose', 0.0) self.cum_na_added += _fresh * enrich.get('Sodium', 0.0) self.cum_cl_added += _fresh * enrich.get('Chloride', 0.0) elif pv > 1e-12: # shrink belly: drain excess dialysate to waste # THE GLUCOSE BOOKING MOVED DOWN into the two branches below (fixed 2026-08-18, # adversarial review). It used to be incremented HERE — four lines before `_to_loop` # decides whether this fluid goes to the SINK or back into the CIRCUIT — so on every # reservoir or PD20_TIDAL_TO_LOOP configuration, glucose that was recirculated and # reused was counted as having gone down the drain, AND counted again downstream. # Measured: the engine's own glucose identity (line ~732) closes to -0.000 mmol on the # default circuit and was broken by 13.5% on every loop-routed one. # This is the SAME defect class already fixed one branch below at "BOOK THE REJECT # AFTER THE STORE, NOT BEFORE" — that fix landed on the reservoir store and missed # this path. Now booked per-transfer, against what actually moved, only when the # destination really is the sink. # WITH A RESERVOIR the surplus goes into the CIRCUIT, so it is regenerated with the # rest of this step's throughput and can be stored clean. Without one it goes # straight to the drain, which is the 14 L/day the attribution audit found. _to_loop = (self.reservoir_cap_L > 0) if self.tidal_to_loop is None \ else bool(self.tidal_to_loop) _dest = self.circuit if _to_loop else self.sink if _fl is None: for L in self.layers: # PD20: proportional withdrawal across the column _mv = (-dV) * (L.volume / pv) # RECORD WHAT MOVED, NOT WHAT WAS ASKED FOR — the same discipline the UF and # lymph legs adopted, missed on this path. transfer_fluid clamps to the source # pool; accumulating the COMMANDED draw over-reported 4.1x on a near-empty # cavity (36 mL held against a 150 mL command). Inert in continuous flow. # Urea booked here too, so the stage attribution covers this waste path. # CONCENTRATIONS READ BEFORE THE TRANSFER (fixed 2026-08-20, adversarial # review of the 2026-08-18 fix). The water counter was corrected to use # transfer_fluid's RETURN, but the two solute counters on the next lines read # L.conc() AFTERWARDS — and Pool.conc() returns 0.0 for an empty pool, so a # CLAMPED transfer (which leaves the layer at exactly 0.0 L) booked the water # and booked the solute as ZERO. That is precisely the near-empty-cavity case # the fix's own comment cites, so the fix addressed half of its own stated # failure. Measured: 1.46 mmol of glucose into the sink booked as 0.000, and # the engine's glucose identity broken by exactly that in ONE step. # The cascade branch below is immune (it leaves a 1e-9 L film), which is why # only this branch was wrong. _cu, _cg = L.conc('Urea'), L.conc('Glucose') _mvd = transfer_fluid(L, _dest, _mv) if not _to_loop: self.cum_dwell_drain_L += _mvd self.cum_urea_dwell_drain += _cu * _mvd self.cum_glu_rejected += _cg * _mvd else: # single bottom lumen: take the LOWEST fluid first, # CLAMP TO WHAT THE CAVITY ACTUALLY HOLDS. The proportional path above cannot # over-drain because it removes a FRACTION (L.volume/pv) of each slab; a cascade # can, because it removes an absolute amount. `dV` is computed from the COMMANDED # volume V_dwell, which drifts from the true volume as ultrafiltration adds # water, so a commanded drain can exceed the real contents. Unclamped this # emptied the cavity to exactly zero and the per-layer transport then divided by # ~0, producing a single-step ultrafiltration excursion of -3160 L. Leave a # residual film so no slab is ever exactly empty. want = min(-dV, max(0.0, pv - 1e-6)) for _Li in self.layers: # cascading up only when a slab is exhausted if want <= 1e-15: break take = min(want, max(0.0, _Li.volume - 1e-9)) if take > 0.0: _tk = transfer_fluid(_Li, _dest, take) # what MOVED, not what was asked if not _to_loop: self.cum_dwell_drain_L += _tk self.cum_urea_dwell_drain += _Li.conc('Urea') * _tk self.cum_glu_rejected += _Li.conc('Glucose') * _tk want -= take self.V_dwell = dwell_L self._apply_dwell_area() if dwell_L is not None: # The EMA must advance on EVERY step the dwell lever is in use, NOT only on steps where # the volume CHANGES. It used to sit inside the branch above: under a constant # non-baseline fill the slew limiter reaches target in ~9 steps, the branch then stops # firing, and the average froze ~8.6% of the way there. control_tbw() subtracts # (V_dwell_avg - fill_L), so the device's own instilled dialysate was permanently # misread as PATIENT fluid — measured 0.91 L of phantom overload at peak fill on a 70 kg # patient, exactly 0 at baseline fill. That biased every constant-fill experiment on the # fill axis, and biased the tidal control arm (rate 0 holds a constant mean, so it froze, # while every cycling arm converged). Found by adversarial review 2026-08-01; pinned by # test_dwell_avg_tracks.py. # PER-MINUTE TIME CONSTANT, NOT PER-STEP (fixed 2026-08-20, adversarial review). # `0.99 * avg + 0.01 * V` has no dt, so its time constant is 100 STEPS — 100 min at # DT=1.0 but 16.7 min at DT=1/6. This estimator feeds control_tbw(), which IS the # regulated variable for both the teacher and the deployed network, so the DEVICE'S OWN # FLUID ESTIMATOR got 5x faster when the integration tick shrank. Measured: 63% step # response 103 min -> 20 min across the swept ticks, and the oscillation the controller # actually sees in control_tbw changed 3.91x (0.232 -> 0.906 L) against the 0.30 L # decision threshold this file names — at the finest tick the estimator noise alone is # three times the threshold. # This is the SAME defect as the §5c glucose one, 630 lines above, and it is the THIRD # tick confound found in the same sweep: patient metabolism, controller gains, and now # the fluid estimator all silently rescaled with dt. # alpha = 1 - 0.99**dt reproduces 0.01 exactly at dt=1 and holds the 100-MINUTE constant # at every other tick, so no DT=1.0 number moves. _alpha = 1.0 - 0.99 ** dt self.V_dwell_avg = (1.0 - _alpha) * self.V_dwell_avg + _alpha * self.V_dwell # Gibbs-Donnan ratio from the IMPERMEANT anion charge (albumin dominates; icodextrin # too where present). r = 1 - |z_protein|*C_protein / (2 * total diffusible ionic charge), # the standard small-charge linearisation. Typical plasma value ~0.96, i.e. serum chloride # sits ~4% BELOW dialysate chloride and serum sodium ~4% above dialysate sodium. _imp = sum(abs(VALENCE.get(_s, 0)) * self.serum.conc(_s, self.serum_vd[_s]) for _s in ('Albumin', 'Globulin', 'Icodextrin')) # denominator is the DIFFUSIBLE SALT concentration (cation charge alone, ~150 mEq/L), # NOT cations+anions — the standard Donnan linearisation r = 1 - z_p*C_p/(2*C_salt). # Using both signs double-counts and gives r ~0.98 where plasma measures ~0.95-0.96. _ion = sum(VALENCE.get(_s, 0) * self.serum.conc(_s, self.serum_vd[_s]) for _s in SOLUTES if VALENCE.get(_s, 0) > 0) _r = 1.0 - _imp / max(1e-9, 2.0 * _ion) _r = min(1.0, max(0.90, _r)) # physiologic bound; never inverts the gradient # ---------- SNAPSHOT (all gradients evaluated at top-of-step) ------- sc = {s: self.serum.conc(s, self.serum_vd[s]) for s in SOLUTES} # (no ICF concentration snapshot: nothing consumes it. When a transcellular flux is added, # read it here as self.tissue.conc(s, self.tissue_vd[s]) — now a TRUE concentration.) lc = [{s: L.conc(s) for s in SOLUTES} for L in self.layers] pv = self._perit_volume() # 1) External IN: oral intake/food (source -> serum) if self.enable['external'] and intake > 0: transfer_fluid(self.source, self.serum, intake * dt, carried_conc={s: intake_c.get(s, 0.0) for s in SOLUTES}) self.cum_na_diet += intake * dt * intake_c.get('Sodium', 0.0) # dietary Na load in self.cum_drink_L += intake * dt # oral water in, L — THIRST-DRIVEN, so it is an # OUTPUT of the simulation, not a prescription # 2) Peritoneal membrane fluxes, per layer, from the snapshot if self.enable['peritoneal']: g0_membrane = self.serum.solutes['Glucose'] # for cum_glu_abs (membrane only) # GEOMETRY: fluid height from the ovoid, then per-slab volume, wall area and depth. # Index 0 is the cavity FLOOR (where the drain lumen sits); index N-1 is the top of the # fluid column (where fresh dialysate enters). # Guard: the ovoid has a finite capacity. If a configuration ever commands a fill # above it, height_at_volume saturates at the top and the slab volumes would silently # be rescaled to match — reporting a geometry that is not the one being simulated. # Fail loudly instead. (Shipped fills are 3.4 L max against 4.0-8.1 L capacity across # the 45-90 kg cohort, so this should never fire; it exists to catch a future change.) # OVERFILL is a physiologic state, not a coding error: the abdominal wall stretches and # intra-abdominal pressure rises steeply. Found 2026-08-07 during pretraining, where a # low-outflow sweep let a 93 kg patient accumulate 8.44 L against an 8.38 L cavity — # the engine previously grew the cavity silently because it had no geometry at all. # _geometry() records the excess; IAP carries it and fatal_check owns the verdict. h_fluid, slabs = self._geometry(pv) # WETTED-VOLUME GUARD (added 2026-08-08). Peritoneal transport must fall to zero as the # cavity empties: with no dialysate there is no exchange. Without this the model treats a # near-empty cavity as having a full membrane facing a ZERO concentration (Pool.conc() # returns 0 below 1e-12 L), so the serum-to-cavity gradient goes to the full serum value # and the fluxes explode. Measured before the guard, simulating a conventional dry-day # cycler that drains the cavity: a single step produced -3160 L of ultrafiltration. # Real peritoneum also retains a residual film after any drain — complete drainage is # not achievable clinically — so scaling smoothly to zero over the last ~150 mL is both # numerically necessary and physiologically right. Inert in all continuous-flow # operation, where the cavity never approaches this volume. wet = min(1.0, max(0.0, pv / PERIT_MIN_WET_L)) self.h_fluid = h_fluid a_tot = sum(sl[2] for sl in slabs) or 1.0 # SNAPSHOT the per-layer volumes that `pv` was measured from. The loop below adds # ultrafiltrate to each layer BEFORE computing that layer's lymph term, so reading # L.volume live while dividing by the pre-loop `pv` mixes a grown numerator with a stale # denominator: sum(vf) then exceeds 1 and total lymphatic uptake exceeds lymph_rate*dt. # Negligible in continuous flow (UF per minute is a few mL against a ~2.45 L cavity, so # sum(vf) is 1 to four decimals) and severe on a near-dry cavity, where a single minute # can ultrafiltrate 475 mL into a 128 mL belly: sum(vf) reached 4.72 and lymph ran 4.6x # its own rate. Snapshotting makes sum(vf) == 1 by construction, at any cavity volume. # (Adversarial review 2026-08-09; the hypotheses of a mid-step refill and of a stale # V_dwell were both tested and refuted before this one.) v_snap = [L.volume for L in self.layers] pv_snap = sum(v_snap) or 1.0 for i, L in enumerate(self.layers): z_mid, v_geo, a_wall = slabs[i] if i < len(slabs) else (0.0, 0.0, 0.0) # HYDROSTATIC: the real standing column above this slab, not a linear stand-in. # 3 cmH2O resting IAP + the water column, converted to mmHg. depth_m = max(0.0, h_fluid - z_mid) # Beyond capacity the wall stretches: 25 cmH2O per litre of excess is the # compliance implied by clinical IAP-vs-fill data (IAP roughly doubles from 2 L to # 3 L of fill in adults). Documented as an order-of-magnitude figure. # INTRAPERITONEAL PRESSURE, recalibrated to measurement 2026-08-09. # Durand PY, Chanliau J, Gamberoni J, Hestin D, Kessler M, Adv Perit Dial 1993;9:46-8 # (PMID 8105960): 34 CAPD patients lying flat at an intraperitoneal volume of # 2820 +/- 319 mL measured IPP 12-14 cmH2O, mean ~13. Durand 1994 (PMID 7999866) # filled 20 adults in 0.5 L steps from 2 to 5 L and found IPP LINEAR in volume, with # 18 cmH2O the maximal acceptable pressure and 20 the stopping rule. # The old law read 7.8 cmH2O at Durand's own 2.82 L — about 40% LOW — which # OVERSTATED the hydrostatic ultrafiltration drive (dP = 17 - IAP) by roughly 51%. # Its slope was 1.1 cmH2O/L across 2-4 L against a measured ~4, and the 25 cmH2O per # litre of OVERFILL was invented: Durand found no steepening at all out to 5 L, so # 25 was ~6x the measured in-range compliance. Both are replaced by the one measured # linear law, anchored so that 2.82 L gives 13 cmH2O. # IAP(V) = IAP_BASE + IAP_SLOPE * V, with the standing water column retained as the # per-slab term that makes the floor slabs see more pressure than the top ones. iap_bulk = IAP_BASE + IAP_SLOPE * (pv + self.overfill_L) iap_i = (iap_bulk + (depth_m * 100.0 - 0.5 * h_fluid * 100.0)) * 0.735 dP = 17.0 - iap_i # Δπ from the ONE shared implementation in physics_core. This used to be a # hand-inlined van't Hoff sum over every solute, a second copy of the osmotic law # that could drift out of step with physics_core.osmotic_pressure_diff — and which # would have done so the moment plasma protein stopped following van't Hoff. # Crystalloids are still van't Hoff; the PROTEIN term is now Landis-Pappenheimer # (see physics_core), which is what closes the transperitoneal oncotic gradient # from 11.2 mmHg to the 22 measured in human PD patients. dpi = osmotic_pressure_diff_conc(lc[i], sc) # AREA-WEIGHTED, not an equal 1/N split: the floor slab is the curved bottom cap # with a high area-to-volume ratio, the upper slabs are wide and shallow-walled. aw = a_wall / a_tot uf = self.Kf * aw * wet * (dP + dpi) * dt # RECORD WHAT MOVED, NOT WHAT WAS ASKED FOR. transfer_fluid clamps to the source # pool and returns the volume it actually moved; accumulating the COMMANDED `uf` and # `lymph_amt` instead reported water that never left the patient. Exact under # continuous flow (nothing ever clamps) but 14.85% high in the cycler comparator, # where a near-dry cavity clamps the lymph and back-filtration legs 80 times a run. # This is also the counter behind the historical -3160 L excursion: the water never # moved, only the accumulator did. (Adversary round 6, S6.) if uf >= 0: # serum -> layer (sieved) uf_moved = transfer_fluid(self.serum, L, uf, carried_conc={s: SIEVING[s] * sc[s] for s in SOLUTES}) else: # back-filtration layer -> serum uf_moved = -transfer_fluid(L, self.serum, -uf, carried_conc=lc[i]) # lymph: layer -> serum (carries layer conc), SCALED BY THE WETTED GUARD like the # ultrafiltration and diffusion legs above. It was not, and because lymph is # apportioned by layer volume FRACTION (sum of vf == 1 by construction), total # lymphatic uptake came out as lymph_rate*dt REGARDLESS of how much fluid the belly # actually held — a cavity down to 100 mL still donated the full 1.07 mL/min, for as # long as you asked. Inert in continuous flow, where the cavity never leaves ~2.45 L # and wet == 1. Badly wrong for a night cycler: through the 15 h dry day the engine # stripped the residual film faster than it could exist, the dwell lever refilled it # from `source` every minute, and that delivered 0.59-0.74 L/day of dialysate nobody # prescribed. That artifact — not physiology — produced the negative net # ultrafiltration in the conventional arm. Found by adversarial review, 2026-08-09. vf = (v_snap[i] / pv_snap) if pv_snap > 1e-12 else 0.0 lymph_amt = self.lymph_rate * vf * wet * dt lymph_moved = transfer_fluid(L, self.serum, lymph_amt, carried_conc=lc[i]) self.cum_uf += uf_moved - lymph_moved # net water off the patient (signed, REALISED) # diffusion serum<->layer (solute only) # GIBBS-DONNAN: albumin cannot cross, so the equilibrium the diffusible ions # relax toward is NOT the raw dialysate concentration. Anions equilibrate to # r x dialysate and cations to dialysate / r, with r < 1 set by the impermeant # protein charge. Uncharged solutes (urea, glucose) are unaffected, r = 1. # PECLET COUPLING (added 2026-08-22, literature audit finding #1). # Diffusion and convection were summed INDEPENDENTLY: the sieved ultrafiltration # above carried solute, and this loop then added the full Fickian flux as if no # solvent were moving. Every published PD kinetic model — Waniewski, and Rippe's # three-pore — attenuates the diffusive term when solvent drag runs the same way, # because the convective sweep flattens the very gradient diffusion feeds on. The # Patlak / Kedem-Katchalsky form is # Js = Jv(1-sigma)*Cbar + PS*(Cb - Cd) * Pe/(e^Pe - 1), # Pe = Jv(1-sigma)/PS # and Pe/(e^Pe - 1) -> 1 as Jv -> 0, so this is EXACTLY the old behaviour at zero # ultrafiltration and only bites when water is actually moving. # Measured before the fix, sodium at PS 4.5 mL/min, S 0.55, Jv 10 mL/min: # Pe = 1.22, factor 0.51 — the engine moved ~2x too much sodium diffusively during # high UF. Urea Pe = 0.34, factor 0.85, ~18% over. That biased sodium removal, the # sodium-dip signature, and EVERY CONDUCTIVITY READING THE CONTROLLER SEES, which # is why this ranked first: it is the controller's primary sensor. # Jv here is the realised water flux for THIS layer this step (L), signed +ve # serum->layer; convert to a rate per unit MTAC to form Pe. Back-filtration # (Jv < 0) opposes diffusion and RAISES the factor above 1, which the same # expression gives without a special case. _jv = (uf_moved - lymph_moved) # net water serum->layer, L this step for s in SOLUTES: z = VALENCE.get(s, 0) tgt = lc[i][s] * (_r if z < 0 else (1.0 / _r if z > 0 else 1.0)) _ps = self.mtac[s] * aw * wet * dt # L this step (permeability-area) if _ps > 1e-15 and abs(_jv) > 1e-15: # Pe USES THE TRANSMISSION COEFFICIENT DIRECTLY, NOT ITS COMPLEMENT. # Kedem-Katchalsky writes Pe = Jv(1-sigma)/PS where sigma is the REFLECTION # coefficient — and in this codebase SIEVING[s] IS (1-sigma): it is used # directly as the convective carriage fraction at :1201 # (`carried_conc = SIEVING[s] * sc[s]`). So (1 - SIEVING) inverts it. # Measured consequence of the inverted form: albumin's SIEVING is 0.01, so # (1-SIEVING) = 0.99 drove its Pe to effectively infinity and switched # protein diffusion OFF. The steady-state transperitoneal oncotic gradient # collapsed from Oberg's measured 22 mmHg to 7.16 — caught by # test_oncotic_landis, which pins exactly that headline claim. _pe = _jv * SIEVING[s] / _ps # Pe/(e^Pe - 1), numerically safe at both tails if abs(_pe) < 1e-8: _f = 1.0 elif _pe > 60.0: _f = 0.0 elif _pe < -60.0: _f = -_pe else: _f = _pe / math.expm1(_pe) else: _f = 1.0 j = _ps * (sc[s] - tgt) * _f transfer_solute(self.serum, L, s, j) # net glucose the membrane moved INTO serum this step (absorption from # the dextrose dialysate): serum gain over the per-layer membrane section. self.cum_glu_abs += self.serum.solutes['Glucose'] - g0_membrane # (REMOVED 2026-08-07) A "WELL-MIXED cavity" homogenizer used to sit here and reset every # layer to the cavity-MEAN concentration on every step. It was the real mixing mechanism # in this engine and it survived the ovoid-geometry rewrite unnoticed, so the vertical # gradient the geometry produced was being ANNIHILATED and re-created from scratch each # minute by that step's inflow/drain alone. Measured: a deliberately imposed 20% gradient # went to 0.000000% in ONE 1-minute step with membrane and flow both off. Any gradient # reported while this block was live was a one-step advection residue, NOT a converged # advection-vs-diffusion balance. Inter-layer mixing is now done ONCE, in the Fickian # block in 4b, where it scales with N through dz = h/N (which is what made the original # fixed-dz Fickian N-dependent and motivated this homogenizer in the first place). # N-independence of clearance is re-verified by test_layer_mixing.py. # 3) Transcellular water PARTITION (serum<->ICF): redistribute water so serum and # ICF are ISO-TONIC at the current effective osmolality — ALGEBRAIC (placed AT # equilibrium each step: no rate constant, cannot oscillate or drain the ICF). # Serum effective osmoles = 2*Na + glucose (impermeant to the cell, σ≈1); ICF # osmoles tracked (self.icf_osmoles). Rising serum glucose raises A -> water # ICF->serum -> translocational (dilutional) hyponatremia (Katz/Hillier). Pure # water (carries NO solute) and serum+tissue volume is preserved, so total body # water and the fluid controller are untouched — only serum CONCENTRATIONS shift. A = 2.0 * self.serum.solutes['Sodium'] + self.serum.solutes['Glucose'] # serum eff osmoles B = self.icf_osmoles # ICF eff osmoles Vpair = self.serum.volume + self.tissue.volume if A + B > 1e-9: Vs_target = A * Vpair / (A + B) shift = Vs_target - self.serum.volume # +ve: tissue->serum (dilute serum) zero = {s: 0.0 for s in SOLUTES} if shift > 1e-12: transfer_fluid(self.tissue, self.serum, shift, carried_conc=zero) elif shift < -1e-12: transfer_fluid(self.serum, self.tissue, -shift, carried_conc=zero) # 4) Apparatus: outflow -> circuit -> REGENERATE -> permeate, enrich, inflow. if self.enable['apparatus'] and q > 0: # DRAIN FROM THE BOTTOM OF THE COLUMN, cascading upward (fixed 2026-08-07). The drain # used to pull the entire q*dt out of layers[0] alone. That silently imposes a # resolution limit: the floor slab is the SMALLEST slab in an ovoid (33 mL of a 2.45 L # fill at N=10) and one minute of drain at 20 mL/min is most of it, so at finer # discretisation the demand exceeds the slab and the transfer clamps. Measured # consequence: BUN 8.23 at N=10 but 10.17 at N=20 — a 23% clearance swing produced # purely by the layer count, which is exactly the N-dependence the (now removed) # well-mixed homogenizer had been papering over. A real drain lumen at the floor # removes the LOWEST fluid in the cavity; if it needs more than the bottom slab holds it # takes the next slab up. Cascading reproduces that and makes clearance grid-independent. glu_before = sum(L.solutes['Glucose'] for L in self.layers) want = q * dt for _Li in self.layers: # bottom -> up: slab 0 is the floor if want <= 1e-15: break take = min(want, max(0.0, _Li.volume - 1e-12)) if take > 0.0: transfer_fluid(_Li, self.circuit, take) want -= take # DO NOT SILENTLY DISCARD AN UNSATISFIABLE DRAIN (added 2026-08-07, adversary finding). # A cavity cannot supply more than it holds, so limiting the drain is correct physics — # but the previous loop dropped the shortfall with no record, so a command of 5 L/min # against a 2.52 L cavity quietly lost 2.48 L of intended drain and every downstream # metric (throughput, cartridge use, clearance) was computed as though it had happened. # Track it so the shortfall is auditable instead of invisible. if want > 1e-12: self.cum_drain_deficit_L += want self.drain_deficit_steps += 1 self.cum_glu_drained += glu_before - sum(L.solutes['Glucose'] for L in self.layers) # glucose that LEFT the belly (not absorbed) # CIRCUIT HOLDUP / transport delay (FULL_PHYSICS_AUDIT §E): the circuit is a # well-mixed device priming/dead volume. Regenerate only the THROUGHPUT this step, # leaving `circuit_holdup` standing, so the regenerated dialysate's composition LAGS # the cavity effluent by ~holdup/q (transit delay). holdup=0 -> feed IS the whole # circuit and it ends empty -> EXACTLY the prior instantaneous path (regression-safe). if self.circuit_holdup > 0.0: feed = Pool(0.0) transfer_fluid(self.circuit, feed, max(0.0, self.circuit.volume - self.circuit_holdup)) else: feed = self.circuit self.cond_out = conductivity(feed) # effluent (drained belly) conductivity # SORBENT UPSTREAM (regen='sorbent+ro'): the column treats the RAW EFFLUENT, before the # reverse-osmosis stage, rather than polishing the permeate after it. The two placements # are NOT equivalent and the difference is capacity, not chemistry: # UPSTREAM sees the whole drained volume at full concentration, so it must be sized # for the entire urea load INCLUDING the fraction the RO would have sent to # waste anyway — capacity spent on solute that was already leaving. # DOWNSTREAM sees only the recovered permeate, and removes exactly the urea that would # otherwise be re-infused into the patient. Same clinical effect per unit of # cartridge, on a much smaller stream. # Upstream does have one genuine advantage: it presents the RO with a cleaner feed, which # matters if the design later runs at higher recovery or drops the recirculation loop. # Measured against each other by exp_sorbent_serial.py. if self.regen == 'sorbent+ro': for _sp, _frac in self.sorbent_strip.items(): _bound = transfer_solute(feed, self.sink, _sp, feed.solutes[_sp] * _frac) self.cum_sorbent_bound[_sp] += _bound if self.regen == 'sorbent': # SORBENT regeneration (REDY / wearable-kidney): strips UREA to waste, recirculates # the rest. NOT USED (device runs RO; sorbent is off the table per standing directive). # UREA-ONLY IS THE DEFAULT AND IT IS THE PROBLEM WITH A CLOSED LOOP. With no reject # stream, urea is the ONLY solute with an exit, so potassium and phosphate build up # behind the bed: measured at 30 days, K 5.39 and phosphate 2.44 mmol/L against # 4.98 and 1.54 on the shipped reverse-osmosis circuit. `sorbent_multi` switches the # recirculating bed to the same PER-SPECIES removal the serial column already uses, # which is what a real multi-layer cartridge does. Default False = prior behaviour. if self.sorbent_multi: for _sp, _frac in self.sorbent_strip.items(): _bound = transfer_solute(feed, self.sink, _sp, feed.solutes[_sp] * _frac) self.cum_sorbent_bound[_sp] += _bound else: for s in ('Urea',): transfer_solute(feed, self.sink, s, feed.solutes[s] * self.sorbent_urea_strip) transfer_fluid(feed, self.permeate, feed.volume) # recirculate all else: # RO split feed -> permeate(clean water) + reject(WASTE incl glucose). osm = 0.0258 * sum(feed.conc(s) * MACRO_OSM_FACTOR.get(s, 1.0) for s in SOLUTES) # CONCENTRATION POLARISATION (film theory) — β=1.0 by DEFAULT, an exact no-op that # reproduces every published number. Set sim.ro_cp_beta > 1 to enable. See # physics_core.cp_beta and exp_polarization.py: at β≥1.1 the 2-pass configuration # stops satisfying the concentration-blind mixer bound, which is why the value matters. beta = cp_beta(osm, self.ro_wp, self.ro_pressure, self.ro_ref_jw, self.ro_cp_beta) jw_raw = self.ro_wp * (self.ro_pressure - beta * osm) jw = max(jw_raw, 0.01) ff = max(0.0, min(1.0, jw_raw / self.ro_ref_jw)) # OSMOTIC WALL (2026-08-13): the commanded recovery is only achievable if the # CONCENTRATE's osmotic pressure stays under the applied pressure. See # ro_max_recovery(). Without this the engine recovered 70% of a feed whose # concentrate would have needed 30.4 bar against 15 applied. r_eff = min(recovery, self.ro_max_recovery(feed, beta)) self.ro_recovery_capped = r_eff < recovery - 1e-12 # OPEN QUESTION, FLAGGED AND DELIBERATELY NOT CHANGED (adversarial review 2026-08-19). # `r_eff` and `ff` are BOTH decreasing functions of the same feed osmolality — r_eff is # the thermodynamic recovery ceiling (1 - pi_feed/P) and ff is that same osmolality # expressed as a flux ratio (jw_raw/ref_jw) — and they are MULTIPLIED here. Measured # live at the shipped point: feed osm 8.655 bar against 15 applied gives r_max 0.4230 # and ff 0.9065, so the engine recovers 0.3835 of the feed, about 9% BELOW the wall it # just solved for. Arguably the two model different limits (thermodynamic ceiling vs # whether the membrane can move that much water in the residence time), which is why # this is not being called a bug. # WHY IT IS NOT BEING "FIXED": the direction is CONSERVATIVE. Under-recovering # OVERSTATES the water bill, so a change here would make the headline number BETTER, # and a reviewer must not quietly improve the number he is reviewing. An adversarial # pass tried to break that reasoning by arguing the defect must FLATTER clearance # (less urea recycled through a membrane that rejects urea at only 0.50 = a wider # gradient) and MEASURED THE OPPOSITE: BUN moves +0.1% at 15 bar and not at all at # the design point. The conservatism holds in both directions that could be measured. # THE SIZE, CORRECTED 2026-08-20 — the earlier note here said "~9% pessimistic ... # 4.31 L/day included" and that attached the wrong number to the wrong operating # point. 9.4% is the recovery loss at the SHIPPED 15-bar circuit; at the 120-bar # design point where 4.31 L/day is measured, ff = (P - pi)/(P - 8) is 0.9942 and the # loss is 0.58%. Measured end to end at the design point the waste effect is 1.5% # (4.469 -> 4.401 L/day), with BUN and absorbed glucose unmoved. The margin is NOT a # stable 9%; it is strongly operating-point dependent, which is itself the reason to # disclose it rather than carry it as a fixed caveat. # NOW DISCLOSED IN THE PUBLISHED SPEC (build_model_spec.py 6.1) rather than living # only in this comment — an adversary's point that "leave it, keep the comment" # quietly converted a live defect in a published document into a note nobody outside # would ever read. It needs the inventor's decision and a re-baseline, not a silent # edit; but the READER is now told. perm_vol = r_eff * ff * feed.volume carried = {s: beta * feed.conc(s) * (1.0 - max(0.0, min(1.0, 1.0 - self.ro_sp[s] / (self.ro_sp[s] + jw)))) for s in SOLUTES} transfer_fluid(feed, self.permeate, perm_vol, carried_conc=carried) self.cum_glu_rejected += feed.solutes['Glucose'] # remaining feed glucose -> waste self.cum_urea_ro1 += feed.solutes['Urea'] # FIRST-pass clearance, booked here self.cum_reject_L += feed.volume transfer_fluid(feed, self.sink, feed.volume) # reject -> waste # RECYCLER (LUK-003 claim 27): re-filter the recovered permeate back through the RO # (ro_passes-1) more times. Each extra pass keeps a cleaner subset of clean water and # rejects the concentrated remainder to waste, so the re-infused dialysate carries less # waste solute (urea) -> wider clearance gradient. ro_passes=1 -> no-op (regression-safe). # SERIAL SORBENT COLUMN on the PERMEATE (regen='ro+sorbent'). Placed here, after # the RO stage and before the recycler, so it polishes the recovered water that is # about to be re-infused. Removal is per-species and fractional per pass. if self.regen == 'ro+sorbent': for _sp, _frac in self.sorbent_strip.items(): _bound = transfer_solute(self.permeate, self.sink, _sp, self.permeate.solutes[_sp] * _frac) self.cum_sorbent_bound[_sp] += _bound for _ in range(max(0, self.ro_passes - 1)): if self.permeate.volume <= 1e-9: break stage = Pool(0.0) transfer_fluid(self.permeate, stage, self.permeate.volume) # permeate -> stage (permeate emptied) osm2 = 0.0258 * sum(stage.conc(s) * MACRO_OSM_FACTOR.get(s, 1.0) for s in SOLUTES) beta2 = cp_beta(osm2, self.ro_wp, self.ro_pressure, self.ro_ref_jw, self.ro_cp_beta) jw2_raw = self.ro_wp * (self.ro_pressure - beta2 * osm2) jw2 = max(jw2_raw, 0.01) ff2 = max(0.0, min(1.0, jw2_raw / self.ro_ref_jw)) # the same wall applies to every subsequent pass, and bites HARDER because the # permeate entering pass 2 has already been concentrated once r2 = min(recovery, self.ro_max_recovery(stage, beta2)) pv2 = r2 * ff2 * stage.volume carried2 = {s: beta2 * stage.conc(s) * (1.0 - max(0.0, min(1.0, 1.0 - self.ro_sp[s] / (self.ro_sp[s] + jw2)))) for s in SOLUTES} transfer_fluid(stage, self.permeate, pv2, carried_conc=carried2) # cleaner subset -> permeate self.cum_glu_rejected += stage.solutes['Glucose'] # rejected glucose -> waste self.cum_urea_ro2p += stage.solutes['Urea'] # SECOND (and later) pass clearance self.cum_reject_L += stage.volume transfer_fluid(stage, self.sink, stage.volume) # this pass reject -> waste self.cond_pre = conductivity(self.permeate) # recovered permeate, PRE-mixer self.cond_in = conductivity(enrich) if enrich else self.cond_pre # infused dialysate # REFILL the cavity back to its physiological dwell volume V_dwell, purchasing fresh # enrichment ONLY for the volume actually infused — never for permeate about to be # discarded. Order (audit fix — was: enrich whole pool THEN dump excess, which counted # the same fresh dextrose in BOTH cum_glu_added and cum_glu_rejected): # (1) dump any EXCESS recovered water UNENRICHED to waste, # (2) enrich the RETAINED permeate up to the target dialysate concentration (deficit- # based = the RECYCLING CREDIT: the RO already recovered solute, fresh only tops the gap), # (3) infuse it, topping up with fresh makeup from source at the target concentration. # Section-2 osmotic UF/lymph already moved water this step, so the cavity sits at # _perit_volume(); net patient fluid removal is exactly the osmotic UF drained as effluent. # PURGE a fixed volume of the regenerated stream, at LOOP concentration, before any of # it is stored or re-infused. This is the loop's glucose (and everything-else) exit. if self.loop_bleed_L_day > 0.0 and self.permeate.volume > 1e-12: _bl = min(self.permeate.volume, self.loop_bleed_L_day * dt / 1440.0) if _bl > 1e-12: self.cum_glu_rejected += self.permeate.solutes['Glucose'] * (_bl / self.permeate.volume) self.cum_urea_bleed += self.permeate.solutes['Urea'] * (_bl / self.permeate.volume) self.cum_loop_bleed_L += _bl transfer_fluid(self.permeate, self.sink, _bl) inflow_target = max(0.0, self.V_dwell - self._perit_volume()) # DRAW ON THE RESERVOIR BEFORE BUYING FRESH. If the regenerator recovered less than the # belly needs — always true under reverse osmosis, where most of the feed is rejected — # stored fluid covers the shortfall first. This is what turns the reservoir from a # buffer into a saving: without it every litre of that shortfall is bought new. if self.reservoir_cap_L > 0 and self.permeate.volume < inflow_target - 1e-12: _pull = min(inflow_target - self.permeate.volume, self.reservoir.volume) if _pull > 1e-12: transfer_fluid(self.reservoir, self.permeate, _pull) self.cum_reservoir_returned_L += _pull if self.permeate.volume > inflow_target + 1e-12: # excess recovered water -> waste (UNENRICHED) excess = self.permeate.volume - inflow_target # STORE WHAT THE BELLY CANNOT TAKE, up to the reservoir's capacity, and discard only # the overflow. The stored fluid is regenerated and unenriched, which is exactly the # state the mixer expects, so returning it later costs only the enrichment deficit. _store = min(excess, max(0.0, self.reservoir_cap_L - self.reservoir.volume)) if _store > 1e-12: transfer_fluid(self.permeate, self.reservoir, _store) self.cum_reservoir_stored_L += _store excess -= _store if excess > 1e-12: # BOOK THE REJECT AFTER THE STORE, NOT BEFORE. This increment used to sit above # the reservoir store and used the PRE-store excess, so glucose parked in the # reservoir — which is returned to the loop and used, never drained — was counted # as having gone to waste. Over-report was exactly permeate_conc * stored_volume, # so it was ZERO at reservoir 0 (today's design point, and the only case # test_glucose_ledger exercises) and grew with reservoir size: every 2 L and # 0.5 L configuration on record over-states its glucose-to-waste. # transfer_fluid moves solute proportionally, so the permeate CONCENTRATION is # unchanged by the store and evaluating it here is exact. self.cum_glu_rejected += self.permeate.solutes['Glucose'] * (excess / self.permeate.volume) self.cum_urea_excess += self.permeate.solutes['Urea'] * (excess / self.permeate.volume) self.cum_excess_dump_L += excess transfer_fluid(self.permeate, self.sink, excess) if enrich and self.permeate.volume > 0: # enrich ONLY the retained permeate for s in SOLUTES: deficit = (enrich.get(s, 0.0) - self.permeate.conc(s)) * self.permeate.volume moved = min(max(0.0, deficit), self.source.solutes[s]) transfer_solute(self.source, self.permeate, s, moved) if s == 'Glucose': self.cum_glu_added += moved # fresh dextrose drawn from supply elif s == 'Sodium': self.cum_na_added += moved # fresh sodium drawn from supply elif s == 'Chloride': self.cum_cl_added += moved # fresh chloride drawn by the mixer take = min(self.permeate.volume, inflow_target) transfer_fluid(self.permeate, self.layers[self.N - 1], take) makeup = inflow_target - take if makeup > 1e-12: # top up with fresh dialysate from source transfer_fluid(self.source, self.layers[self.N - 1], makeup, carried_conc={s: enrich.get(s, 0.0) for s in SOLUTES}) self.cum_glu_added += makeup * enrich.get('Glucose', 0.0) # fresh dextrose in makeup self.cum_na_added += makeup * enrich.get('Sodium', 0.0) # fresh sodium in makeup self.cum_cl_added += makeup * enrich.get('Chloride', 0.0) # fresh chloride in makeup if self.permeate.volume > 1e-12: # numerical remainder -> waste self.cum_glu_rejected += self.permeate.solutes['Glucose'] self.cum_urea_remainder += self.permeate.solutes['Urea'] self.cum_remainder_L += self.permeate.volume transfer_fluid(self.permeate, self.sink, self.permeate.volume) # 4b) Inter-layer VOLUME rebalance. The peritoneal cavity is one connected # fluid: inflow lands on top, outflow drains the bottom, and UF/lymph act # per layer — left alone, the layer volumes diverge (top balloons, bottom # empties so its outflow clamps, and the cavity grows unboundedly). Equalize # layer volumes via paired transfers (carrying solute at source conc), so # the cavity volume is physical and the bottom always has fluid to drain. if self.enable['peritoneal'] or self.enable['apparatus']: # Settle each slab toward its GEOMETRIC volume, not an equal 1/N share. Forcing equal # volumes was itself the mixing mechanism that flattened the cavity: it shuffled fluid # (and with it solute) between layers every step. With ovoid targets the correction is # small — it only removes drift — so a real vertical gradient can persist. pv_now = self._perit_volume() h_now, slabs_t = self._geometry(pv_now) tgt = [sl[1] for sl in slabs_t] if slabs_t else [pv_now / self.N] * self.N ssum = sum(tgt) or 1.0 tgt = [t * pv_now / ssum for t in tgt] # exact conservation of cavity volume # ADJACENT-ONLY settling = physical DISPLACEMENT, not teleportation. The previous # version moved excess from any slab directly to any other, which is instantaneous # perfect mixing: fluid entering at the top appeared throughout the cavity in one step, # and no gradient could survive regardless of the diffusion coefficient (verified — # changing PERIT_MIX_D by 100x moved the gradient by 0.1%). Real fluid entering at the # top DISPLACES the column downward slab by slab. Sweeping adjacent pairs reproduces # that plug-flow displacement while remaining a paired, conservation-safe transfer. # CASCADE form (fixed 2026-08-07). The previous version swept top->down and, whenever a # slab needed MORE than the slab below it held, clamped and STRANDED the deficit — and # because the sweep terminates at i=1, every stranded remainder piled up in slab 0, the # floor. Measured over 2000 steps: worst error 132.7 mL on a 33 mL target slab (399%), # with slabs 2-4 also off by 25-57 mL. Layer volumes were therefore not the geometric # volumes they were claimed to be. # The cascade below is the EXACT plug-flow answer and converges in ONE pass: sweeping # upward, each slab sheds exactly its own surplus (or absorbs its own deficit) across its # top boundary, so slabs 0..N-2 land precisely on target and slab N-1 receives the # remainder, which equals tgt[N-1] because sum(tgt) == pv exactly. It is still # adjacent-only displacement — no slab exchanges with a non-neighbour — and each transfer # carries solute at the donor slab's concentration, as displacement physically does. # ITERATED TO CONVERGENCE (fixed again 2026-08-07, adversary finding). A FIXED two-pass # sweep is not enough. One pass is exact only when no draw is volume-limited; when the # drain has emptied several adjacent slabs, the clamp binds and a deficit refills by only # ONE slab per pass, so two passes cannot repair four empty slabs. Measured on the # two-pass version: at N=40 slab 0 sat at exactly 0.0 mL and at N=60 the bottom FOUR did, # with settle residual 1.5e-2 L. Because Pool.conc() returns 0.0 for an empty pool, those # slabs silently reported ZERO concentration, which flipped the sign of the reported urea # gradient and pinned "spread" at 100%. It was invisible to test_layer_mixing.py because # that test only swept N=4..20, where it does not occur. # Now: sweep until the worst deviation stops improving, bounded by N passes (the most # any deficit can need, since it advances one slab per pass). # MEASURED cost (the earlier claim here that "at N<=30 this exits after the first pass" # was wrong): mean passes/step at nominal flow is 1.000 for N<=20, 1.997 at N=25-30, # 2.993 at N=40, 3.990 at N=60; step time 0.224 ms at N=10 rising to 1.388 ms at N=60. # Published runs use N=10, where it is a single pass, so this is a test-suite cost only. # If the pass budget is ever exhausted the residual is RECORDED rather than dropped — # the same discipline as the drain deficit, and the defect this loop replaced. for _pass in range(max(2, self.N)): worst = 0.0 for i in range(self.N - 1): # bottom -> up, boundary by boundary d_i = self.layers[i].volume - tgt[i] if d_i > 1e-12: # slab too full: push up across boundary transfer_fluid(self.layers[i], self.layers[i + 1], min(d_i, max(0.0, self.layers[i].volume - 1e-9))) elif d_i < -1e-12: # slab short: draw down from above take = min(-d_i, max(0.0, self.layers[i + 1].volume - 1e-9)) if take > 1e-12: transfer_fluid(self.layers[i + 1], self.layers[i], take) worst = max(worst, abs(self.layers[i].volume - tgt[i])) if worst <= 1e-12: break else: _res = max(abs(self.layers[i].volume - tgt[i]) for i in range(self.N)) if _res > 1e-12: self.cum_settle_residual_L += _res self.settle_unconverged_steps += 1 # INTER-LAYER SOLUTE DIFFUSION (added 2026-08-07). Adjacent slabs exchange solute by # Fickian transport across their shared horizontal cross-section. The effective # diffusivity is NOT molecular — the cavity is stirred by respiration, posture and # peristalsis — so PERIT_MIX_D lumps convective mixing with molecular diffusion and is # calibrated to a cavity mixing time of order ten minutes. It is a SENSITIVITY # PARAMETER, not a measured constant: too high and the cavity is one well-mixed pool # (which is what the engine did before this), too low and it stratifies unphysically. # SOLVED IMPLICITLY (backward Euler) ON A TRIDIAGONAL SYSTEM — see the note below on why # neither an explicit flux nor a pairwise relaxation is acceptable here. if (self.enable['peritoneal'] and len(slabs_t) == self.N and h_now > 1e-9 and self.N > 1): dz = h_now / self.N V = [L.volume for L in self.layers] kb = [] # boundary conductance k_{i+1/2}, L/min for i in range(self.N - 1): zb = (i + 1) * dz # shared cross-section at the boundary u = (zb - self.cav[2]) / self.cav[2] a_x = math.pi * self.cav[0] * self.cav[1] * max(0.0, 1.0 - u * u) # m^2 kb.append(PERIT_MIX_D * a_x / max(dz, 1e-9) * 1000.0 * 60.0) if min(V) > 1e-12: # WHY IMPLICIT. The finite-volume balance for the column is # V_i dc_i/dt = k_{i-1/2}(c_{i-1} - c_i) + k_{i+1/2}(c_{i+1} - c_i). # Two earlier forms both failed: # * an EXPLICIT flux k*(ca-cb)*dt, which at these numbers asks for ~13 slab # volumes of exchange per 1-minute step and survived only because an # overshoot cap (measured binding 99.9% of the time) silently replaced the # stated diffusivity as the rate-setter; # * an exact PAIRWISE two-tank relaxation, which is exact for an isolated pair # but is Gauss-Seidel across a chain: its per-pair exponent k*dt/V scales as # N^2, so refining the grid drives every pair to full equilibration and the # cavity collapses to the well-mixed limit. Measured: BUN 8.18 at N=10 but # 10.06 at N=20 — a 23% clearance swing produced purely by the layer count. # Backward Euler on the tridiagonal system is unconditionally stable, has no # cap, converges as the grid is refined, and conserves solute exactly because # its inter-cell fluxes are antisymmetric by construction. Solved with the # Thomas algorithm in O(N); the elimination sweep depends only on the geometry, # so it is done ONCE per step and reused for all 13 species. lo = [0.0] * self.N di = [0.0] * self.N up = [0.0] * self.N for i in range(self.N): kl = kb[i - 1] if i > 0 else 0.0 kr = kb[i] if i < self.N - 1 else 0.0 lo[i] = -dt * kl up[i] = -dt * kr di[i] = V[i] + dt * (kl + kr) cpr = [0.0] * self.N # species-independent forward sweep piv = [0.0] * self.N piv[0] = di[0] cpr[0] = up[0] / piv[0] for i in range(1, self.N): piv[i] = di[i] - lo[i] * cpr[i - 1] cpr[i] = (up[i] / piv[i]) if i < self.N - 1 else 0.0 newc = [0.0] * self.N for sp in SOLUTES: # rhs_i = V_i * c_i^old = the slab's mmol, exactly dpr = [0.0] * self.N dpr[0] = self.layers[0].solutes[sp] / piv[0] for i in range(1, self.N): dpr[i] = (self.layers[i].solutes[sp] - lo[i] * dpr[i - 1]) / piv[i] newc[self.N - 1] = dpr[self.N - 1] for i in range(self.N - 2, -1, -1): newc[i] = dpr[i] - cpr[i] * newc[i + 1] # realise the solution as PAIRED transfers so the closed-universe ledger # stays structurally conservative rather than conservative-by-arithmetic for i in range(self.N - 1): jm = dt * kb[i] * (newc[i] - newc[i + 1]) # +ve: lower -> upper if jm > 0: transfer_solute(self.layers[i], self.layers[i + 1], sp, jm) elif jm < 0: transfer_solute(self.layers[i + 1], self.layers[i], sp, -jm) # 5) External OUT: renal + insensible (serum -> sink) if self.enable['external']: if urine > 0: # PER-SOLUTE renal clearance (was a flat gfr_frac for every solute, which left urine # iso-osmotic to serum and grossly UNDER-cleared urea -> residual-patient BUN biased # high). A concentrating nephron makes urine UREA far above plasma (residual urea # clearance ~ GFR); Na/others stay near the filtered fraction; urine glucose ~0 # (reabsorbed). K is cleared by the distal-secretion term (renal_k) below, so keep # its urine fraction at the filtered level here to avoid double-counting. urine_conc = {s: gfr_frac * sc[s] * URINE_CONC_MULT.get(s, 1.0) for s in SOLUTES} transfer_fluid(self.serum, self.sink, urine * dt, carried_conc=urine_conc) self.cum_patient_loss_L += urine * dt if insensible > 0: transfer_fluid(self.serum, self.sink, insensible * dt, carried_conc={s: 0.0 for s in SOLUTES}) self.cum_patient_loss_L += insensible * dt # 5b) Metabolic urea generation (protein catabolism). Patients continuously # PRODUCE urea; without it BUN->0 and serum osmolality drifts down, tipping # the water/electrolyte balance. Modeled as urea entering serum from the # supply pool (dietary-protein proxy), conservation-correct. ugen = action.get('urea_gen_mmol', 0.0) if ugen > 0: transfer_solute(self.source, self.serum, 'Urea', min(ugen, self.source.solutes['Urea'])) # Creatinine generation (muscle mass): produced at a constant mass-dependent rate, excreted # only by filtration, so it accumulates in anuria until the device removes it — the quantity # D/P creatinine measures. Same source->serum transfer as urea, conservation-correct. cgen = action.get('creatinine_gen_mmol', 0.0) if cgen > 0: transfer_solute(self.source, self.serum, 'Creatinine', min(cgen, self.source.solutes['Creatinine'])) # 5b'') SODIUM diet: FIXED daily dietary salt load (source->serum), diet-driven and # INDEPENDENT of thirst water. Real osmotic thirst drinks ~free water, not saline, so # dietary Na must NOT scale with drinking (the old intake_conc coupling made # hyperglycemia inflate the Na load — see FULL_PHYSICS_AUDIT.md §H). Tracked as dietary Na. dna = action.get('dietary_na_mmol', 0.0) if dna > 0: self.cum_na_diet += transfer_solute(self.source, self.serum, 'Sodium', min(dna, self.source.solutes['Sodium'])) # 5b') POTASSIUM balance: dietary load in (source->serum); fecal/colonic K out # (serum->sink) + stool water out (serum->sink at serum conc). ESRD has no renal K # regulation, so this + the dialysate set serum K. All paired transfers => conserved. dk = action.get('dietary_k_mmol', 0.0) if dk > 0: transfer_solute(self.source, self.serum, 'Potassium', min(dk, self.source.solutes['Potassium'])) stool = action.get('stool_Lmin', 0.0) if stool > 0: # Colon avidly reabsorbs Na/Cl -> fecal fluid Na ~30 mmol/L, NOT serum 140 # (FULL_PHYSICS_AUDIT §J: serum-conc Na gave ~21 mmol/d, too high; real ~5-10). stool_conc = {s: sc[s] for s in SOLUTES} stool_conc['Sodium'] = min(sc['Sodium'], 30.0) stool_conc['Chloride'] = min(sc['Chloride'], 20.0) stool_conc['Albumin'] = 0.0 # stool albumin is zero; the generic serum-conc route stool_conc['Globulin'] = 0.0 # (same for globulin: no protein-losing enteropathy here) # was losing 6.6 g/d of a 66.5 kDa protein transfer_fluid(self.serum, self.sink, stool * dt, carried_conc=stool_conc) self.cum_patient_loss_L += stool * dt fecal_k = action.get('fecal_k_mmol', 0.0) if fecal_k > 0: transfer_solute(self.serum, self.sink, 'Potassium', min(fecal_k, self.serum.solutes['Potassium'])) # Renal DISTAL K SECRETION (serum->sink) — on top of the bulk filtered K in the urine flow # above; net renal K excretion is secretion-dominated (see physics_patient.renal_k_mmol). renal_k = action.get('renal_k_mmol', 0.0) if renal_k > 0: transfer_solute(self.serum, self.sink, 'Potassium', min(renal_k, self.serum.solutes['Potassium'])) # 5c) Glucose metabolism — physiologic two-sided control (conservation-correct): # • HEPATIC GLUCOSE OUTPUT when serum glucose is BELOW the fasting set point: # brisk insulin-INDEPENDENT counter-regulation (glucagon/liver) restores # fasting glucose, prevents the glucose-free-dialysate hypoglycemia trap. # • INSULIN-MEDIATED DISPOSAL when ABOVE set point: uptake scaled by INSULIN # SENSITIVITY Si. Si=1 healthy; T2DM is insulin-RESISTANT (Si~0.2-0.4), so # under the PD dextrose-absorption load serum glucose rises and STAYS high # (the real diabetic hazard). Set point itself is raised for T2DM (fasting # hyperglycemia). Old callers pass glucose_target_mmol -> Si defaults 1.0, # so non-diabetic behavior is unchanged. g_set = action.get('glucose_setpoint_mmol', action.get('glucose_target_mmol', 0.0)) if g_set > 0: si = action.get('insulin_sensitivity', 1.0) want = g_set * self.serum_vd['Glucose'] cur = self.serum.solutes['Glucose'] err = want - cur # RATE, NOT A PER-STEP FRACTION (fixed 2026-08-18, adversarial review). These two terms # carried no `* dt`, so the PATIENT'S OWN liver and insulin ran proportionally faster at a # finer integration tick. Measured on an isolated patient (no device, no controller), a # 25 -> 5 mmol/L relaxation over EXACTLY 60 simulated minutes: 5.92 mmol/L at DT=1.0, # 5.04 at 0.5, 5.000 at 0.25. That is physiology changing with the numerics. # WHY IT MATTERED MORE THAN ITS SIZE: PD20_DT is a live knob and the 2026-08-18 tick sweep # varied it precisely to ask whether a faster CONTROLLER helps. The bias runs in the same # direction as the hypothesis under test, so the fine-tick arms were flattered by a # property of the patient rather than of the controller — the same class of error the # actuator envelopes were added to prevent on the DEVICE side, left unfixed on the # PATIENT side. 0.05 per minute is the rate the DT=1.0 runs have always implied, so this # is an exact no-op at DT=1.0 and every published number is unmoved. # RATE CONSTANT 0.05 -> 0.02 /min, AND HEPATIC OUTPUT IS CAPPED # (fixed 2026-08-22, literature audit finding #6). # 0.05/min is 5%/min. The published intravenous-glucose-tolerance K_G in healthy adults # is 1.5-2.5 %/min, so the model disposed of glucose 2-3x too fast — and with # insulin_sensitivity 0.3 a "T2DM" patient disposed at 1.5 %/min, i.e. at a NORMAL # person's rate. The diabetic arm was therefore not diabetic. Consequence for THIS # device: it suppressed the hyperglycaemia a peritoneal dextrose load actually causes, # and with it the serum-glucose back-pressure that opposes the osmotic gradient — so # ultrafiltration looked easier than it is, exactly where the dextrose lever is worked # hardest. 0.02/min sits at the upper end of the published band. # HEPATIC OUTPUT was also unbounded: err scales with the deficit, so at zero serum # glucose it reached ~3.5 mmol/min against a physiologic endogenous production of # ~2 mg/kg/min = ~0.8 mmol/min for a 70 kg adult. An uncapped liver hides the # hypoglycaemia risk of a glucose-free dialysate, which is a real clinical hazard for # this device. Capped at 2 mg/kg/min, scaled to body mass. _kg = 0.02 if err > 0: # hepatic output (insulin-independent, brisk) _hgo_max = 2.0 / 1000.0 / 180.156 * 1000.0 * self.weight # mmol/min, 2 mg/kg/min prod = min(_kg * err, _hgo_max) * dt transfer_solute(self.source, self.serum, 'Glucose', min(prod, self.source.solutes['Glucose'])) elif err < 0: # insulin-mediated disposal, impaired in T2DM disp = si * _kg * (-err) * dt transfer_solute(self.serum, self.sink, 'Glucose', min(disp, self.serum.solutes['Glucose'])) # 5d) INSULIN-driven transcellular POTASSIUM shift (Na-K-ATPase). Insulin — present # with glucose, scaled by insulin_sensitivity — pumps K serum->ICF; a diabetic (low # Si) shifts LESS and therefore runs a HIGHER serum K. Modeled as relaxation toward an # insulin-set serum/ICF partition: net serum->ICF = gain*(insulin_action*[K]serum − # K_LEAK_REF*[K]icf). At rest (Si=1, glu=5) insulin_action=1 and the two terms cancel # (net 0 — a BUFFER, not a sink). Redistribution only; the K controller resupplies serum. si = action.get('insulin_sensitivity', 1.0) glu_now = self.serum.conc('Glucose', self.serum_vd['Glucose']) insulin_pert = si * 0.08 * max(0.0, glu_now - 5.0) # insulin action ABOVE fasting baseline k_serum = self.serum.conc('Potassium', self.serum_vd['Potassium']) ik = (K_INS_GAIN * insulin_pert * k_serum - K_RETURN * self.ik_displaced) * dt moved = transfer_solute(self.serum, self.tissue, 'Potassium', ik) # +ve serum->ICF self.ik_displaced += moved # cumulative insulin-K displacement (self-returning) self.icf_osmoles += 2.0 * moved # K (+ counter-anion) osmoles follow into/out of ICF # 5f) TRANSCELLULAR POTASSIUM BUFFER (added 2026-08-02). The 3920 mmol intracellular pool # used to be inert: the only serum<->ICF route was the insulin transient above, which # self-returns to zero and cannot answer a serum deficit. Measured consequence: with dietary # K stopped, serum K fell to the fatal floor of 2.0 within ~5 h while ICF K did not move. # Form: a linearised pump-leak. C_eq is a set point that RIDES ON ICF K CONTENT, so at rest # J == 0 exactly (a buffer, never a source or a chronic drain); as the cell pool empties, # C_eq falls and the buffer weakens, which is the physiologic behaviour. # C_eq = C0 + (Q_icf - Q_icf0)/S ; J = K_IC * (C_eq - C_ecf) [+ve = ICF -> serum] # NOTE this is paired with VD_FRACTION['Potassium'] 0.40 -> 0.20: the serum pool used to # carry a whole-body apparent Vd, which WAS the (static) buffer. Changing one without the # other either double-counts or removes buffering entirely. if self.enable['external']: # Same per-solute-Vd basis as the skeletal buffer below. Numerically identical today # (VD_FRACTION['Potassium'] == 0.20 == the ECF reference) but written explicitly so this # buffer cannot silently acquire the phosphate defect if that Vd is ever revised. c_ecf = self.serum.conc('Potassium', self.serum_vd['Potassium']) s_cap = K_BODY_PER_PLASMA_PER_KG * self.weight # weight-scaled capacity k_ic = K_IC_PER_KG * self.weight # weight-scaled speed c_eq = (self._k0_plasma + (self.tissue.solutes['Potassium'] - self._q_icf0) / s_cap) jk = k_ic * (c_eq - c_ecf) * dt # mmol this step, +ve ICF -> serum # NO icf_osmoles adjustment. Adding one here was a bug I introduced on 2026-08-02 and # an adversarial pass caught it: the water partition defines serum effective osmoles as # A = 2*Na_serum + Glu_serum, which does NOT include potassium, so crediting B = # icf_osmoles by 2*mv created osmoles on the ICF side and destroyed none on the serum # side. Vs = A*Vpair/(A+B) then contracted the ECF by up to 36% (5.09 vs 7.92 L in the # 45 kg anuric) and FLIPPED THE SIGN of net sodium retained (-317 vs +87 mmol), while # every reported concentration stayed normal so nothing surfaced it. # Physiologically the right increment is ~zero anyway: transcellular K movement is # largely cation EXCHANGE (Na-K-ATPase 3Na out/2K in; pump-leak efflux charge-balanced # by Na/H entry), so it is close to osmotically silent. The insulin term above keeps its # own +=2.0 increment: that is pre-existing, separately validated behaviour and is NOT # in scope of this fix. if jk > 0: transfer_solute(self.tissue, self.serum, 'Potassium', min(jk, self.tissue.solutes['Potassium'])) elif jk < 0: transfer_solute(self.serum, self.tissue, 'Potassium', min(-jk, self.serum.solutes['Potassium'])) # 5f') SKELETAL BUFFER for calcium, phosphate and magnesium. Bone is a real exchangeable # reservoir on this timescale; without it mineral intake has no sink and serum runs away. if self.enable['external']: for _sp in BONE_CONTENT: # PER-SOLUTE Vd, not raw serum volume. `_min0` stores the initial serum_c values, # which are Vd-based concentrations, so comparing them against solutes/volume mixes # two unit systems. VD_FRACTION['Phosphate'] = 0.30 against the 0.20 ECF reference # made the phosphate comparison run 1.5x high, so the buffer pulled ~230 mmol/day # serum->bone AT REST and reported serum phosphate 15-25% LOW. Ca and Mg were # unaffected only because their VD_FRACTION happens to equal the 0.20 reference — # i.e. it was right by coincidence, not by construction. (Adversary round 6, S3.) _c = self.serum.conc(_sp, self.serum_vd[_sp]) _ceq = self._min0[_sp] + (self.bone.solutes[_sp] - self._bone0[_sp]) / ( BONE_S[_sp] * self.weight / 70.0) _j = BONE_K[_sp] * self.weight / 70.0 * (_ceq - _c) * dt if _j > 0: # bone -> serum transfer_solute(self.bone, self.serum, _sp, min(_j, self.bone.solutes[_sp])) elif _j < 0: # serum -> bone transfer_solute(self.serum, self.bone, _sp, min(-_j, self.serum.solutes[_sp])) # 5f'-bis) PTH-LIKE MINERAL HOMEOSTASIS (added 2026-09-06). The skeletal buffer above is a # FINITE linear buffer: its equilibrium concentration rises as it fills, so serum Ca/Mg/PO4 # drift up with it and never plateau. Real physiology closes that loop with PTH — when serum # Ca/Mg/PO4 rises above its setpoint, PTH falls and the excess is deposited into bone (and # excreted renally where residual function exists). This term is that negative feedback: it # removes the excess into bone at a rate proportional to the excess, so the serum mineral # concentrations stabilise at their setpoints instead of running away. for _sp, _set in (('Calcium', 1.15), ('Magnesium', 0.80), ('Phosphate', 1.6)): _c = self.serum.conc(_sp, self.serum_vd[_sp]) _excess = _c - _set if _excess > 0: _rate = 0.002 * _excess * self.serum_vd[_sp] * dt # ~0.5-day time constant transfer_solute(self.serum, self.sink, _sp, min(_rate, self.serum.solutes[_sp])) # 5f'') NON-OSMOTIC SODIUM STORE (S19/P2). Skin + muscle GAG-bound Na, osmotically INACTIVE: # when serum Na rises the store absorbs Na with ZERO obligated water, buffering tonicity; when # it falls the store releases. Same linearised buffer as bone; zero water carried, so # conservation stays structural through transfer_solute. if self.enable['external'] and SKIN_NA_S > 0: _cna = self.serum.conc('Sodium', self.serum_vd['Sodium']) _ceq = self._na0 + (self.skin.solutes['Sodium'] - self._skin0) / ( SKIN_NA_S * self.weight / 70.0) _kna = self.serum.volume / (SKIN_NA_K_DAYS * 1440.0) # L/min (time constant in days) _j = _kna * (_ceq - _cna) * dt if _j > 0: # skin -> serum (release Na) transfer_solute(self.skin, self.serum, 'Sodium', min(_j, self.skin.solutes['Sodium'])) elif _j < 0: # serum -> skin (deposit Na) transfer_solute(self.serum, self.skin, 'Sodium', min(-_j, self.serum.solutes['Sodium'])) # 5g) NET ENDOGENOUS ACID PRODUCTION (added 2026-08-02). Protein metabolism generates a # fixed acid load that titrates bicarbonate; the engine generated urea from protein but not # the acid. 0.7 mmol/kg/d sits in the sourced 0.7-1.0 band (Remer T, Manz F, PMID 7797810; # Frassetto LA et al., PMID 9734733) and matches dialysis-cohort NEAP of 42.7-58.2 mEq/d # (PMID 26508542) at 70 kg. Buffered acid leaves as exhaled CO2 -> serum HCO3 to sink. neap = action.get('neap_mmol', 0.0) if neap > 0: transfer_solute(self.serum, self.sink, 'Bicarbonate', min(neap, self.serum.solutes['Bicarbonate'])) # ...and the CONJUGATE ANION is RETAINED (added 2026-08-05). Modelling the acid as a # pure bicarbonate sink gave a hyperchloraemic NORMAL-gap acidosis (measured anion gap # 0.5-4.7); ESRD is a HIGH-gap acidosis precisely because sulfate/phosphate/organic # anions accumulate. Without this term Cl + HCO3 must absorb the entire cation charge, # so serum Cl and HCO3 trade off exactly and neither can be brought into range while # the other is — the structural dead end measured in the 2026-08-05 base sweep. transfer_solute(self.source, self.serum, 'Sulfate', min(neap * 0.5, self.source.solutes['Sulfate'])) # divalent: 1 mmol = 2 mEq # 5h) HEPATIC ALBUMIN SYNTHESIS (added 2026-08-02). There was none, so serum albumin fell # 4.6 -> 1.2 g/dL over a 30-day run, its oncotic opposition to ultrafiltration collapsed, # and the controller answered by walking dextrose down — understating published glucose # absorption by ~14%. Restorative first-order (synthesis is up-regulated by # hypoalbuminaemia), capped at the physiologic ~15 g/d maximum. if self.enable['external']: deficit = self._alb0 - self.serum.solutes['Albumin'] if deficit > 0: rate = min(deficit * (0.693 / ALB_SYN_HALFLIFE_MIN), ALB_SYN_MAX_MMOL_MIN) transfer_solute(self.source, self.serum, 'Albumin', min(rate * dt, self.source.solutes['Albumin'])) # 5h-bis) GLOBULIN REPLACEMENT (added 2026-08-09 with the globulin species itself). # Globulin is lost to the effluent like albumin — this engine's transport constants put # that loss at ~2 g/day, which is the right size, since a PD patient's total effluent # protein is 5-15 g/day and albumin is only ~50-60% of it. Without a replacement term # globulin could only ever FALL: an unopposed ~0.6%/day decay strips ~17% of the pool # over a 30-day run and drags the colloid osmotic pressure down with it, which is # exactly the failure mode that made albumin-only synthesis necessary in the first # place (2026-08-02: albumin fell 4.6 -> 1.2 g/dL and the oncotic opposition to # ultrafiltration collapsed). Real PD patients hold globulin at or ABOVE the normal # range indefinitely — plasma cells replace it, and chronic inflammation tends to raise # it — so a restorative term is the physiologic model, not a convenience. # Same first-order restorative form and same time constant as albumin. The ceiling is # scaled by the molecular-weight ratio so the two servos share ONE mass-basis reserve # (3 g/day) rather than the globulin one silently being ~2x as strong in grams. g_deficit = self._glb0 - self.serum.solutes['Globulin'] if g_deficit > 0: g_cap = ALB_SYN_MAX_MMOL_MIN * (PROTEIN_MW_G_PER_MMOL['Albumin'] / PROTEIN_MW_G_PER_MMOL['Globulin']) rate = min(g_deficit * (0.693 / ALB_SYN_HALFLIFE_MIN), g_cap) transfer_solute(self.source, self.serum, 'Globulin', min(rate * dt, self.source.solutes['Globulin'])) # 5i) DIETARY CHLORIDE (added 2026-08-02). There was NONE: 215 mmol/d of cation (Na 165 + # K 50) entered the patient daily with zero accompanying anion, and the DEVICE supplied the # patient's entire chloride. ~90% of dietary sodium is NaCl (IOM, DRI for Water, Potassium, # Sodium, Chloride and Sulfate, 2005, ch. 8); the remaining ~10% is the alkali-precursor # fraction, which pairs with the acid load in 5g. # 5j) DIETARY PHOSPHATE (added 2026-08-03), delivered in MEAL BOLUSES not continuously. # Phosphorus intake tracks protein at ~13-15 mg P per g protein; intestinal absorption is # ~60-70%. Phosphate is also a RETAINED ANION in uraemia and is therefore part of what makes # ESRD a high-anion-gap acidosis — the very species the acid-load model was missing. dphos = action.get('dietary_phos_mmol', 0.0) if dphos > 0: transfer_solute(self.source, self.serum, 'Phosphate', min(dphos, self.source.solutes['Phosphate'])) for _sp, _key in (('Calcium', 'dietary_ca_mmol'), ('Magnesium', 'dietary_mg_mmol')): _amt = action.get(_key, 0.0) if _amt > 0: transfer_solute(self.source, self.serum, _sp, min(_amt, self.source.solutes[_sp])) dcl = action.get('dietary_cl_mmol', 0.0) if dcl > 0: transfer_solute(self.source, self.serum, 'Chloride', min(dcl, self.source.solutes['Chloride'])) # 6) Cori reaction in serum (1 Lactate -> 1 HCO3) if self.enable['cori'] and cori > 0: c = min(cori, self.serum.solutes['Lactate']) self.serum.solutes['Lactate'] -= c self.serum.solutes['Bicarbonate'] += c self.cori_consumed += c self.cori_produced += c # 7) CONSERVATION ASSERTION (closed universe; Cori tracked per species) self._assert_conservation() def _assert_conservation(self, wtol=1e-4, stol=1e-2): """Guard the closed-universe conservation with an ABSOLUTE tolerance. The old relative 1e-7 was multiplied by the total (~1e6 L, dominated by the `source` reservoir) → it only enforced ~0.1 L water / ~14 mmol. Float64 resolves a 1e6-scale sum to ~2e-10 and the real per-run drift is ~1e-8 L (cohort `resid`), so 1e-4 L / 1e-2 mmol is ~1000x tighter yet far above the numerical floor. A hard `raise` (not `assert`) so the guard survives `python -O`.""" # COMPARISONS ARE INVERTED SO NaN RAISES (fixed 2026-08-18, adversarial review). # `abs(nan - x) > tol` is False, so a NaN state sailed through the guard the module docstring # promises will "halt on any violation" — the one corruption that would propagate into every # downstream number. `not (… <= tol)` is True for NaN. The guard was otherwise sound: it # correctly raised on a planted 1 mmol leak. W = total_water(self._pools()) if not (abs(W - self._W0) <= wtol): raise AssertionError(f"WATER non-conservation: {W - self._W0:+.3e} L (tol {wtol})") for s in SOLUTES: tot = total_solute(self._pools(), s) expect = self._S0[s] if s == 'Lactate': expect -= self.cori_consumed elif s == 'Bicarbonate': expect += self.cori_produced if not (abs(tot - expect) <= stol): raise AssertionError(f"{s} non-conservation: {tot - expect:+.3e} mmol (tol {stol})") def _geometry(self, pv): """(h_fluid, slabs) for a cavity volume pv, CACHED on volume to 0.1 mL. slab_geometry runs a composite-Simpson surface integral over every slab and was 57% of step time when called fresh twice per step. Cavity volume moves by microlitres between one-minute steps, so a 0.1 mL cache key is exact for all practical purposes and turns two integrations per step into roughly one per hundred steps. """ # Key on the AXES as well as the volume. Nothing mutates self.cav today, so keying on volume # alone was only latent — but it was proven live-broken by an adversarial probe (doubling # self.cav on a warm cache returned the OLD height unchanged), and a future patient-weight # or cavity-compliance term would silently inherit the previous geometry. Near-free to fix. key = (round(pv, 4), self.cav) if self._geo_cache[0] == key: return self._geo_cache[1], self._geo_cache[2] cap = 4.0 / 3.0 * math.pi * self.cav[0] * self.cav[1] * self.cav[2] * 1000.0 self.overfill_L = max(0.0, pv - cap) h = height_at_volume(min(pv, cap), *self.cav) sl = slab_geometry(h, *self.cav, self.N) self._geo_cache = (key, h, sl) return h, sl def perit_height_m(self): """Standing height of peritoneal fluid above the cavity floor (m), from the ovoid.""" return self._geometry(self._perit_volume())[0] def lumen_pressure(self, z_lumen=None): """Gauge pressure (cmH2O) a transducer in a lumen at height z_lumen would read. THIS IS A DEVICE SIGNAL, not a patient measurement: it senses the standing water column in the cavity and therefore the instilled VOLUME, while touching no blood and reading no concentration. It is the physically available proxy for peritoneal fill. """ z = self.z_outflow if z_lumen is None else z_lumen return lumen_pressure_cmH2O(self.perit_height_m(), z) def volume_from_lumen_pressure(self, p_cmH2O, z_lumen=None, strict=False): """Inverse of lumen_pressure: the volume the device would INFER from a pressure tap. `strict=True` returns None when the fluid line has fallen below the lumen, where the reading is genuinely unresolved (every fill from empty to z_lumen gives the same pressure). It is forwarded rather than swallowed: without it a caller silently receives an UPPER BOUND that can be 13x the true volume for a high lumen on a nearly-empty cavity. The shipped drain lumen sits on the floor (z=0), where the unresolved branch cannot be reached.""" z = self.z_outflow if z_lumen is None else z_lumen return volume_from_pressure_L(p_cmH2O, z, *self.cav, strict=strict) def body_tbw(self): return self.serum.volume + self.tissue.volume + self._perit_volume() def patient_tbw(self): """TRUE patient body water, for ANALYSIS and for the fatal-volume trip — serum + tissue, with the peritoneal cavity EXCLUDED outright. This is the patient's DRAINED weight, which is what a clinician weighs and what "fluid gain" means. Why this exists rather than reusing control_tbw(): control_tbw() is the CONTROLLER's estimate and is deliberately blind — the device may only subtract the dwell volume it BELIEVES it instilled, via V_dwell_avg, a ~100 min EMA. Under continuous flow the fill is constant and the EMA is exact, so the two agree. Under a cycler swinging the belly 0.1 <-> 2.0 L in ~13 min the EMA cannot track, and a dW read at one instant oscillated over a 3.76 L band within a single day against a 0.30 L decision threshold (adversary round 6, S2b) — every conventional- PD verdict computed from it was noise. Excluding the cavity removes the disturbance at the source instead of trying to estimate it. Regression-safe for the continuous-flow arm: dW is always a DIFFERENCE from t=0 and the cavity term there is a constant fill_L, which cancels.""" return self.serum.volume + self.tissue.volume def control_tbw(self): """Patient fluid the weight controller regulates. Same as body_tbw() MINUS the dwell volume the DEVICE has instilled ABOVE the baseline fill — the apparatus knows how much dialysate it put in the belly, so a meal-event dwell BOOST must NOT read as patient fluid overload. At a constant dwell (V_dwell == fill_L) this is IDENTICAL to body_tbw() (regression-safe); during a boost it strips the exogenous excess so the controller doesn't over-ultrafiltrate the patient.""" return self.body_tbw() - (self.V_dwell_avg - self.fill_L) # =========================================================================== # TESTS # =========================================================================== def test_closed_system(): """No externals, no apparatus — only internal peritoneal redistribution. Body TBW + every per-solute total must be invariant to ~1e-9; system drifts toward equilibrium (peritoneal vs serum concentrations converge).""" sim = PhysicsSim(enable={'external': False, 'apparatus': False, 'cori': False}) # Isolate pure SOLUTE redistribution at constant volume: zero the bulk water flows (lymph, UF) # so volumes are fixed and only membrane diffusion plus the inter-layer solve act — the cleanest # internal-conservation + equilibration check. # (CORRECTED 2026-08-07. This comment used to say the equilibration came from "the well-mixed # cavity redistribution" and that "the inter-layer Fickian law is SUPERSEDED and the live engine # never calls it". BOTH halves are now false: the homogenizer was removed, and the live engine # runs a Fickian — a conservative finite-volume diffusion solved implicitly on dz = h/N. The # test passes on membrane diffusion regardless, but the comment was describing an engine that # no longer exists, which is the exact failure test_engine_equivalence2.py exists to prevent.) sim.lymph_rate = 0.0 sim.Kf = 0.0 tbw0 = sim.body_tbw() gap0 = abs(sim.layers[0].conc('Sodium') - sim.serum.conc('Sodium', sim.serum_vd['Sodium'])) for _ in range(2000): sim.step(dt=1.0) # assertion runs inside every step assert abs(sim.body_tbw() - tbw0) < 1e-7, f"body TBW drifted {sim.body_tbw()-tbw0:+.3e}" gap1 = abs(sim.layers[0].conc('Sodium') - sim.serum.conc('Sodium', sim.serum_vd['Sodium'])) assert gap1 < gap0, "diffusion should reduce the serum/peritoneal Na gap" return f"closed system conserved (assertion held 2000 steps); Na gap {gap0:.1f}->{gap1:.2f}" def test_external_balance(): """Only intake + urine. ΔTBW_body must equal intake_in − urine_out exactly.""" sim = PhysicsSim(enable={'peritoneal': False, 'apparatus': False, 'cori': False}) tbw0 = sim.body_tbw() intake_rate = 0.8 / 1440.0 # 0.8 L/day in L/min urine_rate = 0.5 / 1440.0 # 0.5 L/day for _ in range(1440): # 24 h at 1-min steps sim.step(dt=1.0, action={'intake_Lmin': intake_rate, 'intake_conc': {'Sodium': 100.0}, 'urine_Lmin': urine_rate, 'gfr_frac': 0.6}) d_tbw = sim.body_tbw() - tbw0 expected = (intake_rate - urine_rate) * 1440.0 assert abs(d_tbw - expected) < 1e-7, f"ΔTBW {d_tbw:+.4f} != intake-urine {expected:+.4f}" return f"external balance exact: ΔTBW={d_tbw:+.4f} L = intake−urine={expected:+.4f} L (24h)" def test_apparatus_on(): """Full physics, apparatus running. Conservation assertion must hold at EVERY step (it halts otherwise). Report the residual — the legacy code's was 1.75 L/day; this must be floating-point zero.""" sim = PhysicsSim() action = {'outflow_Lmin': 0.02, 'recovery': 0.7, # Cl 97 -> 87, ELECTRONEUTRAL (fixed 2026-08-20). With no K/Ca/Mg in this recipe # the balance is Cl = Na - HCO3 - Lactate = 132 - 35 - 10 = 87; 97 carried a # 10 mEq/L anion excess. The engine's own self-test was specifying a dialysate # it could not infuse, and the actuator-envelope block was silently rewriting it # to 87 on every step — so the test measured one fluid and declared another. # Caught by the new electroneutrality validation the moment it was added, which # is the whole argument for validating instead of correcting. 'enrich': {'Sodium': 132, 'Chloride': 87, 'Bicarbonate': 35, 'Lactate': 10}, 'intake_Lmin': 0.8 / 1440.0, 'intake_conc': {'Sodium': 100.0}, 'urine_Lmin': 0.5 / 1440.0, 'gfr_frac': 0.6, 'insensible_Lmin': 0.7 / 1440.0, 'cori_mmol': 0.001} for _ in range(360): # 6 h at 1-min steps sim.step(dt=1.0, action=action) # explicit residual check over the closed universe. This used to PRINT the residual without # asserting it, so the docstring's "must be floating-point zero" was unenforced and the only # real failure path was step()'s own 1e-4 L guard — a 9e-5 L leak would have passed while the # test reported it. 1e-7 L over 6 h is ~1000x tighter than that guard and still far above the # float64 floor for a ~1e6 L universe (the observed 30-day residual is ~1e-8 L). W = total_water(sim._pools()) resid = W - sim._W0 assert abs(resid) < 1e-7, f"apparatus-on water residual {resid:+.3e} L exceeds 1e-7" return f"apparatus-on 6h: conservation held every step; water residual={resid:+.2e} L (legacy was ~1.75 L/day)" if __name__ == '__main__': tests = [test_closed_system, test_external_balance, test_apparatus_on] ok = True for t in tests: try: print(f" [PASS] {t.__name__}: {t()}") except AssertionError as e: ok = False; print(f" [FAIL] {t.__name__}: {e}") print("RESULT:", "ALL PASS" if ok else "FAILURES") import sys sys.exit(0 if ok else 1)