#!/usr/bin/env python3 """ resistance-evolution.py Science Journaling Club, Volume 2, Issue 1, Fall 2025, "Evolution in Silico". QUESTION -------- Antibiotic resistance is usually argued about as a question of how much drug to use. If you write down a population model with a drug concentration that rises at each dose and decays between doses, the timing and the pattern of dosing turn out to matter as much as the size of the dose. Which schedules actually suppress resistance, and does "finish the course" survive contact with the arithmetic? WHAT THIS PROGRAM IS -------------------- This is a simulation. No bacteria, no patients, no laboratory, no clinical data. The club has computers and people, and the computation IS the experiment. Every number printed below is the output of the model defined in this file. Nothing here is a measurement of a living thing and nothing here is medical advice. MODEL ----- Two strains of one bacterial species share one compartment: S drug-sensitive, maximum net growth psi_max_S R drug-resistant, maximum net growth psi_max_R = psi_max_S * (1 - c) c is the fitness cost of resistance in the absence of drug. Pharmacodynamics follow the sigmoid E_max form of Regoes et al. (2004). For a strain with maximum net growth psi_max, minimum net growth psi_min (< 0) and MIC m, at drug concentration C: psi(C) = psi_max - (psi_max - psi_min) * u / (u - psi_min/psi_max) u = (C/m)^kappa This form is constructed so that psi(0) = psi_max, psi(m) = 0 exactly (that is what MIC means), and psi(inf) = psi_min. Algebraically it collapses to psi(C) = k * psi_max * (1 - u) / (psi_max * u + k), k = -psi_min which is the form used in the code. Pharmacokinetics are one-compartment with instantaneous input and first-order elimination. Each taken dose adds D to C; between doses C decays with rate k_elim = ln2 / t_half. Doses are scheduled every tau hours; under imperfect adherence each scheduled dose is taken independently with probability a. Demography is a stochastic birth-death process, integrated by tau-leaping on a fixed step dt: birth rate b_i = psi_max_i (drug-independent) death rate d_i = (psi_max_i - psi_i(C)) + psi_max_S * N_tot / K so the net growth of strain i is psi_i(C) - psi_max_S * N_tot / K. With no drug and N_tot = K the sensitive strain is exactly at equilibrium and the resistant strain declines at rate c, which is what a fitness cost is. Births are Poisson draws, deaths are binomial draws on the standing population so a population can never go negative. Resistance arises only by mutation: each sensitive birth is resistant with probability mu. VALIDATION PERFORMED BEFORE ANY RESULT IS REPORTED -------------------------------------------------- 1. Deterministic core, stochasticity off. Bisection recovers the concentration at which each strain's net growth crosses zero, and those are compared with the MIC values fed in. The traditional mutant selection window is the closed interval between them. 2. The lower edge of the wider selective window, the minimal selective concentration at which R first out-grows S, is compared against a closed-form solution derived from the pharmacodynamic function (see msc_closed_form). 3. A dense scan of C confirms that R gains on S only above that edge and grows in absolute numbers only below MIC_R. 4. With zero fitness cost the resistant strain never declines with no drug present: the minimum net growth rate over a deterministic drug-free course is printed and must not be negative. 5. With zero mutation supply and no standing resistance, resistance never appears in any stochastic replicate. 6. The mean standing resistant population at the end of the drug-free burn-in is compared with the deterministic mutation-selection balance mu*N/c. 7. The tau-leaping step dt is halved twice on one cell and the outcome probabilities are compared. ASSUMPTIONS, ALL OF THEM STRONG ------------------------------- - One compartment. Real drug reaches different tissues at different concentrations and bacteria do not all sit in plasma. - No immune system. Clearance here is done entirely by the drug. In a real infection the immune response does a great deal of the killing, which makes short courses look better than they do here. - Resistance is one step: one mutation, one fixed MIC shift, one fixed cost. Real resistance is often stepwise, and compensatory mutations erode the cost. - No back mutation, no horizontal gene transfer, no plasmids. - Well mixed. No biofilm, no spatial refuges, no site of infection where drug penetrates poorly. - No persisters, no tolerance, no inoculum effect. - Pharmacokinetics identical in every replicate. Real between-patient variability in clearance is large and would widen every distribution here. - Adherence is independent per dose. Real non-adherence is bursty; people miss runs of doses, not isolated ones. - One patient at a time. Nothing about transmission of resistant strains between hosts, which is the part that matters most for public health. LIMITATIONS THAT BEAR ON THE CONCLUSION --------------------------------------- Because there is no immune clearance, the model asks the drug to do all the work, which biases the comparison toward longer courses. Because resistance here is a single large MIC step, the model cannot show the gradual ratcheting that long sub-inhibitory exposure produces in reality. Both of these are named again in the article. USAGE ----- python resistance-evolution.py > resistance-evolution-output.txt Requires Python 3 and numpy. Master seed is printed at the top of the output and every cell draws an independent stream from it through SeedSequence.spawn, so the whole output is deterministic. """ import math import sys import time import numpy as np # -------------------------------------------------------------------------- # Master seed. Fixed once, never tuned. # -------------------------------------------------------------------------- SEED = 20250915 # -------------------------------------------------------------------------- # Model parameters # -------------------------------------------------------------------------- PSI_MAX_S = 1.0 # per hour, drug-free net growth of the sensitive strain COST = 0.10 # fitness cost of resistance PSI_MAX_R = PSI_MAX_S * (1.0 - COST) PSI_MIN = -4.0 # per hour, maximum kill rate, same for both strains KAPPA = 1.5 # Hill coefficient of the pharmacodynamic curve MIC_S = 1.0 # mg/L MIC_R = 16.0 # mg/L K_CAP = 1.0e9 # carrying capacity, cells T_HALF = 3.0 # hours, elimination half-life K_ELIM = math.log(2.0) / T_HALF MU = 1.0e-9 # resistant mutants per sensitive birth DT = 0.02 # hours, tau-leaping step during treatment DT_BURN = 0.05 # hours, step during the drug-free burn-in BURN_H = 48.0 # hours of drug-free burn-in before treatment starts FOLLOW_H = 168.0 # hours of drug-free observation after the last dose FIX_FRAC = 0.9 # N_R above this fraction of K with S gone is absorbing FIX_MIN = 1.0e6 # cells, floor for calling resistance established OUT = sys.stdout def line(ch="-", n=78): print(ch * n) def head(title): print() line("=") print(title) line("=") # -------------------------------------------------------------------------- # Pharmacodynamics # -------------------------------------------------------------------------- def psi(C, psi_max, mic, kappa=KAPPA, psi_min=PSI_MIN): """Net growth rate at concentration C. psi(0)=psi_max, psi(mic)=0.""" k = -psi_min u = np.power(np.asarray(C, dtype=float) / mic, kappa) return k * psi_max * (1.0 - u) / (psi_max * u + k) def msc_closed_form(p_s, p_r, mic_s, mic_r, kappa=KAPPA, psi_min=PSI_MIN): """Concentration where psi_R(C) = psi_S(C), solved in closed form. Writing u_S = (C/A)^kappa = alpha X, u_R = (C/B)^kappa = beta X with X = C^kappa, alpha = A^-kappa, beta = B^-kappa, and k = -psi_min, setting k p (1 - alpha X) / (p alpha X + k) = k q (1 - beta X) / (q beta X + k) and clearing denominators gives a linear equation in X whose solution is X = k (p - q) / (p q alpha - p q beta + k p alpha - k q beta) with C = X^(1/kappa). p = psi_max_S, q = psi_max_R. """ k = -psi_min a = mic_s ** (-kappa) b = mic_r ** (-kappa) denom = p_s * p_r * a - p_s * p_r * b + k * p_s * a - k * p_r * b X = k * (p_s - p_r) / denom return X ** (1.0 / kappa) def bisect(f, lo, hi, tol=1e-14, itmax=400): flo, fhi = f(lo), f(hi) if flo * fhi > 0: raise ValueError("bisection bracket does not straddle a root") for _ in range(itmax): mid = 0.5 * (lo + hi) fm = f(mid) if flo * fm <= 0: hi, fhi = mid, fm else: lo, flo = mid, fm if hi - lo < tol: break return 0.5 * (lo + hi) # -------------------------------------------------------------------------- # Deterministic core: the same rate laws with the randomness switched off # -------------------------------------------------------------------------- def det_rates(ns, nr, C, psi_max_r=PSI_MAX_R, mu=MU): ntot = ns + nr dens = PSI_MAX_S * ntot / K_CAP gs = psi(C, PSI_MAX_S, MIC_S) - dens gr = psi(C, psi_max_r, MIC_R) - dens dns = gs * ns - mu * PSI_MAX_S * ns dnr = gr * nr + mu * PSI_MAX_S * ns return dns, dnr def det_run(ns0, nr0, conc_fn, hours, dt=0.01, psi_max_r=PSI_MAX_R, mu=MU): """RK4 on the deterministic limit. Returns time series.""" n = int(round(hours / dt)) ts = np.empty(n + 1) ns_a = np.empty(n + 1) nr_a = np.empty(n + 1) ns, nr = float(ns0), float(nr0) ts[0], ns_a[0], nr_a[0] = 0.0, ns, nr for i in range(n): t = i * dt def f(tt, a, b): return det_rates(a, b, conc_fn(tt), psi_max_r, mu) k1 = f(t, ns, nr) k2 = f(t + dt / 2, ns + dt / 2 * k1[0], nr + dt / 2 * k1[1]) k3 = f(t + dt / 2, ns + dt / 2 * k2[0], nr + dt / 2 * k2[1]) k4 = f(t + dt, ns + dt * k3[0], nr + dt * k3[1]) ns += dt / 6 * (k1[0] + 2 * k2[0] + 2 * k3[0] + k4[0]) nr += dt / 6 * (k1[1] + 2 * k2[1] + 2 * k3[1] + k4[1]) ns = max(ns, 0.0) nr = max(nr, 0.0) ts[i + 1], ns_a[i + 1], nr_a[i + 1] = t + dt, ns, nr return ts, ns_a, nr_a # -------------------------------------------------------------------------- # Stochastic engine # -------------------------------------------------------------------------- def run_cell(seed_seq, n_rep, dose, tau_h, n_doses, adherence, follow_h=FOLLOW_H, mu=MU, psi_max_r=PSI_MAX_R, cost_free=False, dt=DT, burn_h=BURN_H, record=False, nr0=0, k_cap=K_CAP, mic_r=MIC_R, immune=0.0): """Simulate n_rep independent treatment courses. Returns a dict of per-replicate outcomes plus, optionally, a recorded trajectory for replicate 0. Outcome codes: 0 cleared, 1 resistance established, 2 still infected. """ rng = np.random.default_rng(seed_seq) pmr = PSI_MAX_S if cost_free else psi_max_r NS = np.full(n_rep, int(k_cap), dtype=np.int64) NR = np.full(n_rep, int(nr0), dtype=np.int64) C = np.zeros(n_rep, dtype=np.float64) idx = np.arange(n_rep) outcome = np.full(n_rep, 2, dtype=np.int8) t_clear = np.full(n_rep, np.nan) max_NR = np.zeros(n_rep, dtype=np.int64) NR_at_start = np.zeros(n_rep, dtype=np.int64) NS_end = np.zeros(n_rep, dtype=np.int64) NR_end = np.zeros(n_rep, dtype=np.int64) rec_t, rec_C, rec_S, rec_R = [], [], [], [] def step(dtl, decay): nonlocal NS, NR, C ntot = NS + NR dens = PSI_MAX_S * ntot / k_cap uS = np.power(C / MIC_S, KAPPA) uR = np.power(C / mic_r, KAPPA) k = -PSI_MIN psiS = k * PSI_MAX_S * (1.0 - uS) / (PSI_MAX_S * uS + k) psiR = k * pmr * (1.0 - uR) / (pmr * uR + k) dS = (PSI_MAX_S - psiS) + dens + immune dR = (pmr - psiR) + dens + immune bS = rng.poisson(PSI_MAX_S * NS * dtl) bR = rng.poisson(pmr * NR * dtl) if mu > 0.0: mut = rng.binomial(bS, mu) else: mut = np.zeros_like(bS) xS = rng.binomial(NS, -np.expm1(-dS * dtl)) xR = rng.binomial(NR, -np.expm1(-dR * dtl)) NS = NS + bS - mut - xS NR = NR + bR + mut - xR C = C * decay # ---- burn-in, no drug ------------------------------------------------- nburn = int(round(burn_h / DT_BURN)) for _ in range(nburn): step(DT_BURN, 1.0) NR_at_start[idx] = NR max_NR[idx] = np.maximum(max_NR[idx], NR) # ---- treatment and follow-up ----------------------------------------- steps_per_dose = int(round(tau_h / dt)) total_h = n_doses * tau_h + follow_h nsteps = int(round(total_h / dt)) decay = math.exp(-K_ELIM * dt) if adherence >= 1.0: taken = np.ones((n_rep, n_doses), dtype=bool) else: taken = rng.random((n_rep, n_doses)) < adherence rec_every = max(1, int(round(0.1 / dt))) # record every 0.1 h for s in range(nsteps): t = s * dt j, rem = divmod(s, steps_per_dose) if rem == 0 and j < n_doses: C = C + dose * taken[idx, j] if record and (s % rec_every == 0) and len(idx) > 0 and idx[0] == 0: rec_t.append(t) rec_C.append(float(C[0])) rec_S.append(int(NS[0])) rec_R.append(int(NR[0])) step(dt, decay) max_NR[idx] = np.maximum(max_NR[idx], NR) ext = (NS + NR) == 0 if ext.any(): new = ext & np.isnan(t_clear[idx]) if new.any(): t_clear[idx[new]] = t + dt if (s % 25) == 24 and not record: # recorded single-replicate runs are never compacted away, so the # trajectory keeps being written for the whole observation window fixed = (NS == 0) & (NR >= FIX_FRAC * k_cap) done = ext | fixed if done.any(): orig = idx[done] outcome[orig] = np.where(ext[done], 0, 1) NS_end[orig] = NS[done] NR_end[orig] = NR[done] keep = ~done NS, NR, C, idx = NS[keep], NR[keep], C[keep], idx[keep] if len(idx) == 0: break if len(idx) > 0: NS_end[idx] = NS NR_end[idx] = NR cleared = (NS + NR) == 0 resist = (~cleared) & (NR >= min(FIX_MIN, 0.01 * k_cap)) & (NR > NS) o = np.full(len(idx), 2, dtype=np.int8) o[resist] = 1 o[cleared] = 0 outcome[idx] = o return dict(outcome=outcome, t_clear=t_clear, max_NR=max_NR, NR_at_start=NR_at_start, NS_end=NS_end, NR_end=NR_end, rec=(np.array(rec_t), np.array(rec_C), np.array(rec_S), np.array(rec_R))) def summarise(res, n_rep): o = res["outcome"] p_ok = float((o == 0).mean()) p_rs = float((o == 1).mean()) p_pr = float((o == 2).mean()) se = lambda p: math.sqrt(max(p * (1 - p), 0.0) / n_rep) tc = res["t_clear"][np.isfinite(res["t_clear"])] return dict(p_success=p_ok, se_success=se(p_ok), p_resist=p_rs, se_resist=se(p_rs), p_persist=p_pr, se_persist=se(p_pr), n_clear=int(len(tc)), t_clear_mean=float(tc.mean()) if len(tc) else float("nan"), t_clear_sd=float(tc.std(ddof=1)) if len(tc) > 1 else float("nan"), t_clear_med=float(np.median(tc)) if len(tc) else float("nan")) def conc_schedule(dose, tau_h, n_doses): def f(t): if t < 0: return 0.0 n = min(int(t // tau_h) + 1, n_doses) tot = 0.0 for j in range(n): tot += dose * math.exp(-K_ELIM * (t - j * tau_h)) return tot return f # -------------------------------------------------------------------------- # Small statistics helpers, written out rather than imported # -------------------------------------------------------------------------- def wilson(k, n, z=1.959964): """95% Wilson score interval for a binomial proportion. Used alongside the plain standard error because several cells of the grid come back at exactly 0 or exactly 1, where the normal interval has zero width and is simply wrong. """ if n == 0: return float("nan"), float("nan") p = k / n d = 1.0 + z * z / n c = (p + z * z / (2.0 * n)) / d h = z * math.sqrt(p * (1.0 - p) / n + z * z / (4.0 * n * n)) / d return max(0.0, c - h), min(1.0, c + h) def rankdata(x): x = np.asarray(x, dtype=float) order = np.argsort(x, kind="mergesort") r = np.empty(len(x), dtype=float) r[order] = np.arange(1, len(x) + 1, dtype=float) vals, inv, cnt = np.unique(x, return_inverse=True, return_counts=True) for i, c in enumerate(cnt): if c > 1: m = inv == i r[m] = r[m].mean() return r def spearman(a, b): ra, rb = rankdata(a), rankdata(b) if ra.std() == 0 or rb.std() == 0: return float("nan") return float(np.corrcoef(ra, rb)[0, 1]) def pk_indices(dose, tau_h, msc, mic_s=MIC_S, mic_r=MIC_R): """Steady-state pharmacokinetic summaries for a regimen, drug only.""" r = math.exp(-K_ELIM * tau_h) cpk = dose / (1.0 - r) grid = np.linspace(0.0, tau_h, 20001) ct = cpk * np.exp(-K_ELIM * grid) return dict(cmax=cpk, cmin=cpk * r, fT_s=float((ct > mic_s).mean()), fT_r=float((ct > mic_r).mean()), fT_win=float(((ct > msc) & (ct < mic_r)).mean()), auc_tau=dose / K_ELIM) # ========================================================================== # MAIN # ========================================================================== def report(tag, sm, n): ks = int(round(sm["p_success"] * n)) kr = int(round(sm["p_resist"] * n)) lo_s, hi_s = wilson(ks, n) lo_r, hi_r = wilson(kr, n) return ("%s %7.4f +-%.4f [%.4f,%.4f] %7.4f +-%.4f [%.4f,%.4f]" % (tag, sm["p_success"], sm["se_success"], lo_s, hi_s, sm["p_resist"], sm["se_resist"], lo_r, hi_r)) def main(): t_start = time.time() print("THE DOSING SCHEDULE DECIDES WHETHER RESISTANCE WINS") print("Science Journaling Club, Volume 2 Issue 1, Fall 2025") print("Everything below is simulation output. No bacteria were involved,") print("no patients were involved, and none of it is medical advice.") line("=") print("master seed : %d" % SEED) print("python : %s" % sys.version.split()[0]) print("numpy : %s" % np.__version__) line("-") print("psi_max_S %8.4f /h psi_max_R %8.4f /h (cost c = %.3f)" % (PSI_MAX_S, PSI_MAX_R, COST)) print("psi_min %8.4f /h kappa %8.4f" % (PSI_MIN, KAPPA)) print("MIC_S %8.4f mg/L MIC_R %8.4f mg/L (ratio %.1f)" % (MIC_S, MIC_R, MIC_R / MIC_S)) print("K %8.3g cells mu %8.3g /birth" % (K_CAP, MU)) print("t_half %8.4f h k_elim %8.6f /h" % (T_HALF, K_ELIM)) print("dt %8.4f h burn-in %8.1f h follow-up %6.1f h" % (DT, BURN_H, FOLLOW_H)) ss = np.random.SeedSequence(SEED) streams = iter(ss.spawn(500)) # ---------------------------------------------------------------------- head("VALIDATION 1 MUTANT SELECTION WINDOW, DETERMINISTIC CORE") # ---------------------------------------------------------------------- print("Stochasticity off. Bisection to 1e-14 on the net growth rate of each") print("strain at density zero. The concentration where a strain's net growth") print("crosses zero is that strain's MIC by definition, so recovering the two") print("input MICs is a check that the pharmacodynamic function was coded") print("correctly and that the window edges are where we think they are.") print() f_s = lambda c: float(psi(c, PSI_MAX_S, MIC_S)) f_r = lambda c: float(psi(c, PSI_MAX_R, MIC_R)) mic_s_found = bisect(f_s, 1e-6, 1000.0) mic_r_found = bisect(f_r, 1e-6, 1000.0) print("%-34s %16s %16s %14s" % ("quantity", "club value", "input/analytic", "difference")) line("-") print("%-34s %16.12f %16.12f %14.2e" % ("lower MSW edge, psi_S(C)=0", mic_s_found, MIC_S, mic_s_found - MIC_S)) print("%-34s %16.12f %16.12f %14.2e" % ("upper MSW edge, psi_R(C)=0", mic_r_found, MIC_R, mic_r_found - MIC_R)) print() print("Mutant selection window produced by the model:") print(" [%.12f, %.12f] mg/L, width %.9f mg/L" % (mic_s_found, mic_r_found, mic_r_found - mic_s_found)) print("Mutant selection window fed in:") print(" [%.12f, %.12f] mg/L, width %.9f mg/L" % (MIC_S, MIC_R, MIC_R - MIC_S)) print("agreement to %.1e mg/L, which is the bisection tolerance." % max(abs(mic_s_found - MIC_S), abs(mic_r_found - MIC_R), 1e-15)) # ---------------------------------------------------------------------- head("VALIDATION 2 MINIMAL SELECTIVE CONCENTRATION, CLOSED FORM") # ---------------------------------------------------------------------- print("The traditional window starts at MIC_S because that is where the") print("sensitive strain stops growing. Relative selection for resistance") print("starts earlier, at the concentration where psi_R first exceeds psi_S.") print("The closed form is derived in msc_closed_form() in this file and was") print("solved by hand before the code was written.") print() msc_bis = bisect(lambda c: f_r(c) - f_s(c), 1e-9, MIC_R) msc_cf = msc_closed_form(PSI_MAX_S, PSI_MAX_R, MIC_S, MIC_R) print("%-34s %16s %16s %14s" % ("quantity", "club value", "input/analytic", "difference")) line("-") print("%-34s %16.12f %16.12f %14.2e" % ("MSC, bisection vs closed form", msc_bis, msc_cf, msc_bis - msc_cf)) print("%-34s %16.12f %16.12f %14.2e" % ("MSC in units of MIC_S", msc_bis / MIC_S, msc_cf / MIC_S, (msc_bis - msc_cf) / MIC_S)) print() print("Selective window produced by the model : [%.6f, %.6f] mg/L" % (msc_bis, mic_r_found)) print("that is : [MIC_S/%.2f, %.0f x MIC_S]" % (MIC_S / msc_bis, mic_r_found / MIC_S)) print("Gullberg et al. (2011) measured minimal selective concentrations of") print("MIC/4 to MIC/230 for real drug-strain pairs. Ours sits at MIC/%.1f," % (MIC_S / msc_bis)) print("inside that published range but at its weak-selection end.") # ---------------------------------------------------------------------- head("VALIDATION 3 DENSE SCAN OF THE WINDOW") # ---------------------------------------------------------------------- print("Scanning C from 0 to 40 mg/L in 40001 steps and asking two questions") print("at each point: does R gain on S (psi_R > psi_S), and does R actually") print("grow in absolute numbers (psi_R > 0)? Both must be true inside the") print("window and at no point outside it.") print() cs = np.linspace(0.0, 40.0, 40001) pS = psi(cs, PSI_MAX_S, MIC_S) pR = psi(cs, PSI_MAX_R, MIC_R) gains = pR > pS grows = pR > 0.0 inside = (cs > msc_cf) & (cs < mic_r_found) bad_gain = int(np.sum(gains != (cs > msc_cf))) bad_grow = int(np.sum(grows != (cs < mic_r_found))) bad_both = int(np.sum((gains & grows) != inside)) print("points scanned : %d" % len(cs)) print("points where R gains on S but C <= MSC : %d" % int(np.sum(gains & (cs <= msc_cf)))) print("points where R does not gain but C > MSC : %d" % int(np.sum((~gains) & (cs > msc_cf)))) print("mismatches, 'R gains' vs C > MSC : %d" % bad_gain) print("mismatches, 'R grows' vs C < MIC_R : %d" % bad_grow) print("mismatches, 'R gains and grows' vs in window : %d" % bad_both) print("verdict : %s" % ("PASS" if (bad_gain == 0 and bad_grow == 0 and bad_both == 0) else "FAIL")) print() print("%10s %14s %14s %14s %12s" % ("C (mg/L)", "psi_S (/h)", "psi_R (/h)", "psi_R-psi_S", "in window")) line("-") for c in [0.0, 0.05, 0.1, msc_cf, 0.5, 1.0, 2.0, 4.0, 8.0, 16.0, 24.0, 40.0]: a = float(psi(c, PSI_MAX_S, MIC_S)) b = float(psi(c, PSI_MAX_R, MIC_R)) inw = "yes" if (msc_cf < c < mic_r_found) else "no" print("%10.6f %14.6f %14.6f %14.6f %12s" % (c, a, b, b - a, inw)) # ---------------------------------------------------------------------- head("VALIDATION 4 ZERO FITNESS COST, DRUG ABSENT") # ---------------------------------------------------------------------- print("With c = 0 the resistant strain must never decline when there is no") print("drug. Deterministic run, 500 h, C = 0, N_S(0) = 1e8, N_R(0) = 1e3,") print("cost set to zero, mutation switched off.") print() ts0, ns0, nr0 = det_run(1e8, 1e3, lambda t: 0.0, 500.0, dt=0.02, psi_max_r=PSI_MAX_S, mu=0.0) dnr = np.diff(nr0) print("%-38s : %s" % ("minimum change in N_R over any step", "%.6e cells" % float(dnr.min()))) print("%-38s : %s" % ("minimum net growth rate of R", "%.6e /h" % float((dnr / 0.02 / np.maximum(nr0[:-1], 1e-30)).min()))) print("%-38s : %.6e -> %.6e" % ("N_R at t = 0 / t = 500 h", nr0[0], nr0[-1])) print("%-38s : %s" % ("verdict", "PASS, R never declines" if dnr.min() >= -1e-9 else "FAIL")) print() print("For contrast, the same run with the real cost c = %.2f:" % COST) ts1, ns1, nr1 = det_run(1e8, 1e3, lambda t: 0.0, 500.0, dt=0.02, psi_max_r=PSI_MAX_R, mu=0.0) print("%-38s : %.6e -> %.6e" % ("N_R at t = 0 / t = 500 h", nr1[0], nr1[-1])) print("R declines once the population fills the niche, which is what a cost") print("is. Decay rate over the last 100 h: %.6f /h against the input %.6f." % (-math.log(nr1[-1] / nr1[-5001]) / 100.0, COST)) # ---------------------------------------------------------------------- head("VALIDATION 5 ZERO MUTATION SUPPLY") # ---------------------------------------------------------------------- print("mu = 0, N_R(0) = 0, 2000 stochastic replicates on a schedule that") print("breeds resistance readily when mutation is on (D = 8 mg/L, tau = 24 h,") print("7 days). One resistant cell anywhere means the bookkeeping is wrong.") print() n_v5 = 2000 r5 = run_cell(next(streams), n_v5, dose=8.0, tau_h=24.0, n_doses=7, adherence=1.0, mu=0.0, nr0=0) print("%-42s : %d" % ("replicates", n_v5)) print("%-42s : %d" % ("largest N_R seen in any replicate", int(r5["max_NR"].max()))) print("%-42s : %d" % ("replicates with any resistant cell ever", int((r5["max_NR"] > 0).sum()))) print("%-42s : %d" % ("replicates classified as resistance", int((r5["outcome"] == 1).sum()))) print("%-42s : %.4f" % ("P(treatment success) with no mutation", float((r5["outcome"] == 0).mean()))) print("%-42s : %s" % ("verdict", "PASS" if r5["max_NR"].max() == 0 else "FAIL")) # ---------------------------------------------------------------------- head("VALIDATION 6 STANDING RESISTANCE AT MUTATION-SELECTION BALANCE") # ---------------------------------------------------------------------- print("Before treatment the model runs %g h of drug-free burn-in from" % BURN_H) print("N_S = K, N_R = 0. The deterministic balance is mu*K/c, approached with") print("relaxation time 1/c = %g h, so the accepted value at %g h is the" % (1.0 / COST, BURN_H)) print("finite-time solution mu*K/c * (1 - exp(-c t)).") print() bal_inf = MU * K_CAP / COST bal = bal_inf * (1.0 - math.exp(-COST * BURN_H)) r6 = run_cell(next(streams), 4000, dose=0.0, tau_h=1.0, n_doses=1, adherence=1.0, follow_h=0.0) nrs = r6["NR_at_start"].astype(float) m = nrs.mean() sem = nrs.std(ddof=1) / math.sqrt(len(nrs)) z = (m - bal) / sem print("%-38s %14s %16s %12s" % ("quantity", "club value", "input/analytic", "difference")) line("-") print("%-38s %14.4f %16.4f %12.4f" % ("mean standing N_R at dose 1", m, bal, m - bal)) print("%-38s %14s %16.4f %12s" % ("asymptotic balance mu*K/c", "-", bal_inf, "-")) print("%-38s : %.4f cells" % ("standard error of the mean", sem)) print("%-38s : %.2f" % ("difference in standard errors", z)) print("%-38s : %d of %d (%.2f%%)" % ("replicates with zero standing R", int((nrs == 0).sum()), len(nrs), 100.0 * (nrs == 0).mean())) print("%-38s : %.1f cells" % ("median standing N_R", float(np.median(nrs)))) print("%-38s : %.1f cells" % ("variance of standing N_R", float(nrs.var(ddof=1)))) print("%-38s : %.1f cells" % ("analytic variance lambda*d/(d-b)^2", MU * K_CAP * PSI_MAX_S / (COST ** 2))) if abs(z) < 3.0: print("%-38s : PASS, within 3 standard errors" % "verdict") else: print("%-38s : DISAGREEMENT, %.2f standard errors out" % ("verdict", z)) # ---------------------------------------------------------------------- head("VALIDATION 7 TAU-LEAPING STEP SIZE") # ---------------------------------------------------------------------- print("One cell (D = 8 mg/L, tau = 12 h, 7 days, full adherence) run at three") print("step sizes with 2000 replicates each. A leap that is too coarse makes") print("the outcome probabilities drift with dt.") print() print("%8s %8s %16s %16s %10s" % ("dt (h)", "n", "P(success)", "P(resistance)", "seconds")) line("-") dt_rows = [] for d in (0.08, 0.04, 0.02): t0 = time.time() rr = run_cell(next(streams), 2000, dose=8.0, tau_h=12.0, n_doses=14, adherence=1.0, dt=d) sm = summarise(rr, 2000) el = time.time() - t0 dt_rows.append((d, sm)) print("%8.3f %8d %9.4f+-%.4f %9.4f+-%.4f %10.1f" % (d, 2000, sm["p_success"], sm["se_success"], sm["p_resist"], sm["se_resist"], el)) da, db = dt_rows[-2][1], dt_rows[-1][1] dz = (da["p_success"] - db["p_success"]) / math.sqrt( da["se_success"] ** 2 + db["se_success"] ** 2) print() print("P(success) at dt = 0.04 against dt = 0.02 differs by %.4f, which is" % (da["p_success"] - db["p_success"])) print("%.2f standard errors. dt = %.2f h is used for everything below." % (dz, DT)) # ---------------------------------------------------------------------- head("EXPERIMENT 1 DOSE SIZE BY DOSING INTERVAL, PERFECT ADHERENCE") # ---------------------------------------------------------------------- N_E1 = 2000 doses = [2.0, 4.0, 8.0, 16.0, 32.0] taus = [4.0, 6.0, 8.0, 12.0, 24.0] course_h = 168.0 print("7-day course, adherence 1.0, %d replicates per cell, %d cells," % (N_E1, len(doses) * len(taus))) print("%d simulated courses. D is the concentration jump per dose in mg/L," % (N_E1 * len(doses) * len(taus))) print("so with MIC_S = 1 mg/L the column is also D/MIC_S. Brackets are 95%") print("Wilson score intervals, which stay honest when a cell reads 0 or 1.") print() print("%5s %5s %5s %8s %-40s %-40s %8s" % ("D", "tau", "n", "AUC", "P(success) +-SE [95% Wilson]", "P(resistance) +-SE [95% Wilson]", "t_clear")) line("-") e1 = {} tw = {} for D in doses: for tau in taus: nd = int(round(course_h / tau)) rr = run_cell(next(streams), N_E1, dose=D, tau_h=tau, n_doses=nd, adherence=1.0) sm = summarise(rr, N_E1) ix = pk_indices(D, tau, msc_bis) auc = nd * D / K_ELIM e1[(D, tau)] = (sm, auc, nd, rr) tw[(D, tau)] = ix ks = int(round(sm["p_success"] * N_E1)) kr = int(round(sm["p_resist"] * N_E1)) ls, hs = wilson(ks, N_E1) lr, hr = wilson(kr, N_E1) tc = sm["t_clear_med"] print("%5.1f %5.1f %5d %8.1f %7.4f +-%.4f [%.4f,%.4f] " "%7.4f +-%.4f [%.4f,%.4f] %8s" % (D, tau, nd, auc, sm["p_success"], sm["se_success"], ls, hs, sm["p_resist"], sm["se_resist"], lr, hr, ("%.2f h" % tc) if np.isfinite(tc) else "-")) print() print("Pharmacokinetic summaries of the same 25 regimens at steady state,") print("from the concentration curve alone with no population model in them:") print("%5s %5s %9s %9s %10s %10s %12s" % ("D", "tau", "Cmax", "Cmin", "fT>MIC_S", "fT>MIC_R", "fT in window")) line("-") for D in doses: for tau in taus: ix = tw[(D, tau)] print("%5.1f %5.1f %9.3f %9.3f %10.4f %10.4f %12.4f" % (D, tau, ix["cmax"], ix["cmin"], ix["fT_s"], ix["fT_r"], ix["fT_win"])) print() print("Which pharmacokinetic index predicts the outcome? Spearman rank") print("correlation of each index against P(success) across all 25 cells.") print() keys = sorted(e1.keys()) psucc = np.array([e1[k][0]["p_success"] for k in keys]) presist = np.array([e1[k][0]["p_resist"] for k in keys]) cand = { "total AUC over the course": np.array([e1[k][1] for k in keys]), "AUC per interval / MIC_S": np.array([tw[k]["auc_tau"] / MIC_S for k in keys]), "Cmax / MIC_S": np.array([tw[k]["cmax"] / MIC_S for k in keys]), "Cmax / MIC_R": np.array([tw[k]["cmax"] / MIC_R for k in keys]), "fraction of time above MIC_S": np.array([tw[k]["fT_s"] for k in keys]), "fraction of time above MIC_R": np.array([tw[k]["fT_r"] for k in keys]), "fraction of time in the window": np.array([tw[k]["fT_win"] for k in keys]), } print("%-34s %14s %14s" % ("index", "rho vs P(ok)", "rho vs P(res)")) line("-") ranked = [] for nm, v in cand.items(): a = spearman(v, psucc) b = spearman(v, presist) ranked.append((abs(a), nm, a, b)) print("%-34s %14.4f %14.4f" % (nm, a, b)) ranked.sort(reverse=True) print() print("Best single predictor of treatment success: %s (rho = %+.4f)." % (ranked[0][1], ranked[0][2])) print("Total drug delivered over the course has rho = %+.4f, so on this grid" % spearman(cand["total AUC over the course"], psucc)) print("how much drug you give predicts the outcome far less well than when") print("the concentration sits above the resistant strain's MIC.") print() print("Mechanism check. Splitting every replicate in the grid by whether any") print("resistant cell was standing at the moment the first dose was taken:") print() tot0 = tot0ok = totp = totpok = 0 for k in keys: rr = e1[k][3] z0 = rr["NR_at_start"] == 0 tot0 += int(z0.sum()) tot0ok += int((rr["outcome"][z0] == 0).sum()) totp += int((~z0).sum()) totpok += int((rr["outcome"][~z0] == 0).sum()) print("replicates starting with zero standing resistance : %d" % tot0) print(" of which treatment succeeded : %d (%.4f)" % (tot0ok, tot0ok / max(tot0, 1))) print("replicates starting with resistance already there : %d" % totp) print(" of which treatment succeeded : %d (%.4f)" % (totpok, totpok / max(totp, 1))) print("odds ratio : %.2f" % ((tot0ok / max(tot0 - tot0ok, 1)) / max(totpok / max(totp - totpok, 1), 1e-12))) # ---------------------------------------------------------------------- head("EXPERIMENT 2 MATCHED TOTAL EXPOSURE") # ---------------------------------------------------------------------- N_E2 = 2000 print("Dosing interval fixed at 8 h. The product (number of doses) x (dose") print("size) is held at %.0f mg/L, so every row delivers the same total drug" % 288.0) print("and the same area under the concentration curve. Only the shape of the") print("schedule changes. %d replicates per row." % N_E2) print() PROD = 288.0 plans = [2, 3, 4, 5, 7, 10, 14] print("%6s %7s %9s %9s %-40s %-40s" % ("days", "doses", "D", "AUC", "P(success) +-SE [95% Wilson]", "P(resistance) +-SE [95% Wilson]")) line("-") e2 = {} for days in plans: nd = int(round(days * 24.0 / 8.0)) D = PROD / nd rr = run_cell(next(streams), N_E2, dose=D, tau_h=8.0, n_doses=nd, adherence=1.0) sm = summarise(rr, N_E2) auc = nd * D / K_ELIM e2[days] = (sm, D, nd, auc) ks, kr = int(round(sm["p_success"] * N_E2)), int(round(sm["p_resist"] * N_E2)) ls, hs = wilson(ks, N_E2) lr, hr = wilson(kr, N_E2) print("%6d %7d %9.3f %9.1f %7.4f +-%.4f [%.4f,%.4f] " "%7.4f +-%.4f [%.4f,%.4f]" % (days, nd, D, auc, sm["p_success"], sm["se_success"], ls, hs, sm["p_resist"], sm["se_resist"], lr, hr)) print() best = max(e2.items(), key=lambda kv: kv[1][0]["p_success"]) worst = min(e2.items(), key=lambda kv: kv[1][0]["p_success"]) d1, d2 = best[1][0], worst[1][0] zz = (d1["p_success"] - d2["p_success"]) / math.sqrt( d1["se_success"] ** 2 + d2["se_success"] ** 2) print("Best row : %2d days at D = %7.3f, P(success) = %.4f" % (best[0], best[1][1], d1["p_success"])) print("Worst row : %2d days at D = %7.3f, P(success) = %.4f" % (worst[0], worst[1][1], d2["p_success"])) print("Difference: %.4f, which is %.1f standard errors, on identical drug." % (d1["p_success"] - d2["p_success"], zz)) print("Resistance runs from %.4f (%d days) to %.4f (%d days), a factor of %.0f." % (min(v[0]["p_resist"] for v in e2.values()), best[0], max(v[0]["p_resist"] for v in e2.values()), worst[0], max(v[0]["p_resist"] for v in e2.values()) / max(min(v[0]["p_resist"] for v in e2.values()), 1e-9))) # ---------------------------------------------------------------------- head("EXPERIMENT 3 COURSE LENGTH AT A FIXED DOSE") # ---------------------------------------------------------------------- N_E3 = 2000 print("The matched-exposure ladder changes two things at once, the dose and") print("the length. This experiment changes only the length. D = 32 mg/L every") print("8 h, stopped after 1, 2, 3, 5, 7, 10 or 14 days. Total drug therefore") print("rises fourteen-fold down the table. %d replicates per row." % N_E3) print() print("%6s %7s %9s %-40s %-40s %10s" % ("days", "doses", "AUC", "P(success) +-SE [95% Wilson]", "P(resistance) +-SE [95% Wilson]", "med t_clear")) line("-") e3 = {} for days in [1, 2, 3, 5, 7, 10, 14]: nd = int(round(days * 24.0 / 8.0)) rr = run_cell(next(streams), N_E3, dose=32.0, tau_h=8.0, n_doses=nd, adherence=1.0) sm = summarise(rr, N_E3) auc = nd * 32.0 / K_ELIM e3[days] = (sm, nd, auc) ks, kr = int(round(sm["p_success"] * N_E3)), int(round(sm["p_resist"] * N_E3)) ls, hs = wilson(ks, N_E3) lr, hr = wilson(kr, N_E3) print("%6d %7d %9.1f %7.4f +-%.4f [%.4f,%.4f] " "%7.4f +-%.4f [%.4f,%.4f] %10.2f h" % (days, nd, auc, sm["p_success"], sm["se_success"], ls, hs, sm["p_resist"], sm["se_resist"], lr, hr, sm["t_clear_med"])) print() a1, a14 = e3[1][0], e3[14][0] zz = (a14["p_success"] - a1["p_success"]) / math.sqrt( a1["se_success"] ** 2 + a14["se_success"] ** 2) print("1 day against 14 days on the same dose: P(success) %.4f vs %.4f," % (a1["p_success"], a14["p_success"])) print("a difference of %.4f, which is %.2f standard errors. The extra 13 days" % (a14["p_success"] - a1["p_success"], zz)) print("cost %.0f mg/L of exposure, %.0f times the first day's." % (e3[14][2] - e3[1][2], (e3[14][2] - e3[1][2]) / e3[1][2])) tcs = [e3[d][0]["t_clear_med"] for d in e3] print("Median time to clearance across all rows: %.2f to %.2f h. In this" % (min(tcs), max(tcs))) print("model the drug finishes inside the first day and every later dose is") print("acting on a population that is already gone.") # ---------------------------------------------------------------------- head("EXPERIMENT 4 MISSED DOSES") # ---------------------------------------------------------------------- N_E4 = 2000 print("Adherence is only interesting on a regimen that works when taken") print("properly, so the reference is D = 32 mg/L every 8 h, the schedule that") print("cleared every replicate in Experiment 1. Each scheduled dose is taken") print("independently with probability a. %d replicates per level, two course" % N_E4) print("lengths, 3 days and 7 days.") print() e4 = {} for days in (3, 7): nd = int(round(days * 24.0 / 8.0)) print("course length %d days (%d scheduled doses)" % (days, nd)) print("%7s %9s %10s %-40s %-40s" % ("a", "missed", "AUC taken", "P(success) +-SE [95% Wilson]", "P(resistance) +-SE [95% Wilson]")) line("-") for a in (1.0, 0.95, 0.9, 0.85, 0.8, 0.7, 0.6, 0.5): rr = run_cell(next(streams), N_E4, dose=32.0, tau_h=8.0, n_doses=nd, adherence=a) sm = summarise(rr, N_E4) e4[(days, a)] = sm ks = int(round(sm["p_success"] * N_E4)) kr = int(round(sm["p_resist"] * N_E4)) ls, hs = wilson(ks, N_E4) lr, hr = wilson(kr, N_E4) print("%7.2f %9.2f %10.1f %7.4f +-%.4f [%.4f,%.4f] " "%7.4f +-%.4f [%.4f,%.4f]" % (a, nd * (1 - a), nd * a * 32.0 / K_ELIM, sm["p_success"], sm["se_success"], ls, hs, sm["p_resist"], sm["se_resist"], lr, hr)) print() s10, s08 = e4[(7, 1.0)], e4[(7, 0.8)] zz = (s10["p_success"] - s08["p_success"]) / math.sqrt( s10["se_success"] ** 2 + s08["se_success"] ** 2 + 1e-24) print("Seven-day course, adherence 1.00 -> 0.80: P(success) changes by %.4f" % (s08["p_success"] - s10["p_success"])) print("(%.2f standard errors) and P(resistance) by %+.4f." % (zz, s08["p_resist"] - s10["p_resist"])) t10, t08 = e4[(3, 1.0)], e4[(3, 0.8)] print("Three-day course, adherence 1.00 -> 0.80: P(success) changes by %.4f" % (t08["p_success"] - t10["p_success"])) print("and P(resistance) by %+.4f." % (t08["p_resist"] - t10["p_resist"])) print() print("Control. The same average amount of drug delivered on a schedule with") print("no gaps in it: every dose taken, but each dose scaled down by a. If") print("adherence hurt only by cutting the total, these two would match.") print("%7s %-34s %-34s" % ("a", "ragged schedule (missed doses)", "smooth schedule (smaller doses)")) line("-") for a in (0.8, 0.7, 0.6, 0.5): nd = 21 rr = run_cell(next(streams), N_E4, dose=32.0 * a, tau_h=8.0, n_doses=nd, adherence=1.0) sm = summarise(rr, N_E4) rg = e4[(7, a)] print("%7.2f ok %.4f res %.4f ok %.4f res %.4f" % (a, rg["p_success"], rg["p_resist"], sm["p_success"], sm["p_resist"])) e4[("smooth", a)] = sm # ---------------------------------------------------------------------- head("ARITHMETIC WHAT ONE MISSED DOSE BUYS THE RESISTANT STRAIN") # ---------------------------------------------------------------------- print("No simulation in this section. It is the pharmacokinetics and the") print("pharmacodynamic function integrated by hand, to show where the missed") print("dose results come from. Elimination rate k = ln2/t_half = %.6f /h." % K_ELIM) print() def t_fall(c_from, c_to): return math.log(c_from / c_to) / K_ELIM def log_fold(psi_max, mic, c0, hours, n=400000): tt = np.linspace(0.0, hours, n + 1) cc = c0 * np.exp(-K_ELIM * tt) return float(np.trapezoid(psi(cc, psi_max, mic), tt)) for tag, D0, TAU0 in (("short sharp arm", 32.0, 8.0), ("long mild arm", 288.0 / 42.0, 8.0)): ix0 = pk_indices(D0, TAU0, msc_bis) cmx = ix0["cmax"] print("%s: D = %.4f mg/L every %.0f h" % (tag, D0, TAU0)) print(" steady-state peak Cmax = D / (1 - exp(-k tau)) = %.4f mg/L" % cmx) print(" steady-state trough Cmin = %.4f mg/L" % ix0["cmin"]) for nm, lvl in (("MIC_R", mic_r_found), ("MIC_S", mic_s_found), ("MSC", msc_bis)): if cmx > lvl: tf = t_fall(cmx, lvl) print(" time from peak down to %-5s = %8.4f h (%.4f of a %.0f h interval)" % (nm, tf, min(tf / TAU0, 1.0), TAU0)) else: print(" time from peak down to %-5s = never above it" % nm) for hrs, lab in ((TAU0, "on schedule"), (2 * TAU0, "one dose missed")): lr = log_fold(PSI_MAX_R, MIC_R, cmx, hrs) ls = log_fold(PSI_MAX_S, MIC_S, cmx, hrs) print(" %-16s (%4.1f h): resistant x %12.4e, sensitive x %12.4e" % (lab, hrs, math.exp(lr), math.exp(ls))) lr1 = log_fold(PSI_MAX_R, MIC_R, cmx, TAU0) lr2 = log_fold(PSI_MAX_R, MIC_R, cmx, 2 * TAU0) print(" the gap multiplies the resistant strain by a further %.4e" % math.exp(lr2 - lr1)) if lr2 > 0: need = math.log(1.0e6 / max(bal, 1e-12)) / lr2 print(" gaps needed to take %.1f standing mutants to 1e6 cells: %.2f" % (bal, need)) print() print("Reading those two lines together: on the short sharp schedule a") print("resistant cell is worse off at the end of every interval than at the") print("start, and a single missed dose reverses that. On the long mild") print("schedule it gains ground on every interval whether or not a dose is") print("missed, which is why missing doses barely changes that arm at all.") # ---------------------------------------------------------------------- head("EXPERIMENT 5 MONTE CARLO CONVERGENCE") # ---------------------------------------------------------------------- N_E5 = 20000 print("One cell run long, D = 8 mg/L every 12 h for 7 days, full adherence,") print("%d replicates, with the running estimate recorded as trials pile up." % N_E5) print() t0 = time.time() r5c = run_cell(next(streams), N_E5, dose=8.0, tau_h=12.0, n_doses=14, adherence=1.0) el5 = time.time() - t0 sm5 = summarise(r5c, N_E5) o5 = r5c["outcome"] print("P(success) = %.5f +- %.5f" % (sm5["p_success"], sm5["se_success"])) print("P(resistance) = %.5f +- %.5f" % (sm5["p_resist"], sm5["se_resist"])) print("P(persist) = %.5f +- %.5f" % (sm5["p_persist"], sm5["se_persist"])) print("wall clock = %.1f s" % el5) print() print("%8s %12s %12s %12s %12s %12s" % ("n", "P(success)", "se", "P(resist)", "se", "|p_n - p_N|/se")) line("-") marks = [100, 250, 500, 1000, 2000, 4000, 8000, 12000, 16000, 20000] conv_rows = [] pN = sm5["p_success"] for n in marks: ps = float((o5[:n] == 0).mean()) pr = float((o5[:n] == 1).mean()) es = math.sqrt(max(ps * (1 - ps), 1e-12) / n) er = math.sqrt(max(pr * (1 - pr), 1e-12) / n) conv_rows.append((n, ps, es, pr, er)) print("%8d %12.5f %12.5f %12.5f %12.5f %12.2f" % (n, ps, es, pr, er, abs(ps - pN) / es)) print() rse0 = 100 * conv_rows[0][2] / max(conv_rows[0][1], 1e-12) rseN = 100 * conv_rows[-1][2] / max(conv_rows[-1][1], 1e-12) print("Relative standard error on P(success) falls from %.2f%% at n = 100" % rse0) print("to %.2f%% at n = %d, a factor of %.2f for a factor of %.0f in trials." % (rseN, N_E5, rse0 / rseN, N_E5 / 100)) print("The 1/sqrt(n) prediction for that factor is %.2f." % math.sqrt(N_E5 / 100)) # ---------------------------------------------------------------------- head("EXPERIMENT 6 WHERE A DIFFERENT MODELLING CHOICE CHANGES THE ANSWER") # ---------------------------------------------------------------------- N_E6 = 1000 print("Two arms carried through every sensitivity run: the short sharp arm") print("(D = 32 mg/L q8h for 3 days) and the long mild arm (D = 6.857 mg/L") print("q8h for 14 days). They deliver identical total drug. %d replicates" % N_E6) print("per cell. The baseline row is the model as specified above.") print() print("%-42s %-22s %-22s" % ("variant", "short: ok / resistance", "long: ok / resistance")) line("-") variants = [ ("baseline", {}), ("fitness cost c = 0.02", dict(psi_max_r=PSI_MAX_S * 0.98)), ("fitness cost c = 0.30", dict(psi_max_r=PSI_MAX_S * 0.70)), ("fitness cost c = 0.50", dict(psi_max_r=PSI_MAX_S * 0.50)), ("MIC_R = 4 x MIC_S", dict(mic_r=4.0)), ("MIC_R = 8 x MIC_S", dict(mic_r=8.0)), ("MIC_R = 32 x MIC_S", dict(mic_r=32.0)), ("inoculum K = 1e7", dict(k_cap=1.0e7)), ("inoculum K = 1e8", dict(k_cap=1.0e8)), ("mutation rate mu = 1e-10", dict(mu=1.0e-10)), ("mutation rate mu = 1e-8", dict(mu=1.0e-8)), ("immune clearance 0.10 /h", dict(immune=0.10)), ("immune clearance 0.30 /h", dict(immune=0.30)), ("immune clearance 0.50 /h", dict(immune=0.50)), ] e6 = {} for nm, kw in variants: ra = run_cell(next(streams), N_E6, dose=32.0, tau_h=8.0, n_doses=9, adherence=1.0, **kw) rb = run_cell(next(streams), N_E6, dose=288.0 / 42.0, tau_h=8.0, n_doses=42, adherence=1.0, **kw) sa, sb = summarise(ra, N_E6), summarise(rb, N_E6) e6[nm] = (sa, sb) print("%-42s %8.4f / %8.4f %8.4f / %8.4f" % (nm, sa["p_success"], sa["p_resist"], sb["p_success"], sb["p_resist"])) print() base = e6["baseline"] print("Baseline gap in P(success), short minus long : %+.4f" % (base[0]["p_success"] - base[1]["p_success"])) gaps = [(nm, v[0]["p_success"] - v[1]["p_success"]) for nm, v in e6.items()] gaps.sort(key=lambda t: t[1]) print("Smallest gap across all variants : %+.4f (%s)" % (gaps[0][1], gaps[0][0])) print("Largest gap across all variants : %+.4f (%s)" % (gaps[-1][1], gaps[-1][0])) flips = [nm for nm, g in gaps if g <= 0.0] print("Variants where the short arm stops winning : %s" % (", ".join(flips) if flips else "none")) # ---------------------------------------------------------------------- head("EXPERIMENT 7 TWO TRAJECTORIES RECORDED IN FULL") # ---------------------------------------------------------------------- print("Single replicates recorded every 0.1 h so the article can draw them.") print("Panel A: D = 32 mg/L every 8 h, 3 days. Peaks above MIC_R.") print("Panel B: D = 8 mg/L every 24 h, 7 days. Sits inside the window.") print() trA = run_cell(next(streams), 1, dose=32.0, tau_h=8.0, n_doses=9, adherence=1.0, record=True, follow_h=48.0) trB = run_cell(next(streams), 1, dose=8.0, tau_h=24.0, n_doses=7, adherence=1.0, record=True, follow_h=48.0) for nm, tr in (("A", trA), ("B", trB)): t, c, s_, r_ = tr["rec"] o = int(tr["outcome"][0]) lab = {0: "cleared", 1: "resistance established", 2: "still infected"}[o] print("panel %s: standing N_R at dose 1 = %d, outcome = %s" % (nm, int(tr["NR_at_start"][0]), lab)) print("panel %s: max N_R = %d, final N_S = %d, final N_R = %d" % (nm, int(tr["max_NR"][0]), int(tr["NS_end"][0]), int(tr["NR_end"][0]))) print("panel %s: peak C = %.3f mg/L, last recorded C = %.3f mg/L" % (nm, float(c.max()), float(c[-1]))) if np.isfinite(tr["t_clear"][0]): print("panel %s: cleared at t = %.2f h" % (nm, float(tr["t_clear"][0]))) print() print("DATA BLOCK trajectoryA t_hours C_mgL N_S N_R (every 0.5 h)") tA, cA, sA, rA = trA["rec"] stride = 5 for i in range(0, len(tA), stride): print(" %7.2f %10.4f %14d %14d" % (tA[i], cA[i], sA[i], rA[i])) print("DATA BLOCK trajectoryB t_hours C_mgL N_S N_R (every 0.5 h)") tB, cB, sB, rB = trB["rec"] for i in range(0, len(tB), stride): print(" %7.2f %10.4f %14d %14d" % (tB[i], cB[i], sB[i], rB[i])) # ---------------------------------------------------------------------- head("MACHINE-READABLE BLOCKS FOR THE FIGURES AND THE LATEX") # ---------------------------------------------------------------------- print("BLOCK grid D tau P_ok SE P_res SE P_persist SE fT_s fT_r fT_win Cmax") for D in doses: for tau in taus: sm, auc, nd, _ = e1[(D, tau)] ix = tw[(D, tau)] print(" %.1f %.1f %.5f %.5f %.5f %.5f %.5f %.5f %.4f %.4f %.4f %.3f" % (D, tau, sm["p_success"], sm["se_success"], sm["p_resist"], sm["se_resist"], sm["p_persist"], sm["se_persist"], ix["fT_s"], ix["fT_r"], ix["fT_win"], ix["cmax"])) print("BLOCK matched days D n_doses AUC P_ok SE P_res SE") for days in sorted(e2): sm, D, nd, auc = e2[days] print(" %d %.4f %d %.2f %.5f %.5f %.5f %.5f" % (days, D, nd, auc, sm["p_success"], sm["se_success"], sm["p_resist"], sm["se_resist"])) print("BLOCK duration days n_doses AUC P_ok SE P_res SE") for days in sorted(e3): sm, nd, auc = e3[days] print(" %d %d %.2f %.5f %.5f %.5f %.5f" % (days, nd, auc, sm["p_success"], sm["se_success"], sm["p_resist"], sm["se_resist"])) print("BLOCK adherence days a P_ok SE P_res SE") for days in (3, 7): for a in (1.0, 0.95, 0.9, 0.85, 0.8, 0.7, 0.6, 0.5): sm = e4[(days, a)] print(" %d %.2f %.5f %.5f %.5f %.5f" % (days, a, sm["p_success"], sm["se_success"], sm["p_resist"], sm["se_resist"])) print("BLOCK smooth a P_ok SE P_res SE") for a in (0.8, 0.7, 0.6, 0.5): sm = e4[("smooth", a)] print(" %.2f %.5f %.5f %.5f %.5f" % (a, sm["p_success"], sm["se_success"], sm["p_resist"], sm["se_resist"])) print("BLOCK convergence n P_ok SE P_res SE") for n, ps, es, pr, er in conv_rows: print(" %d %.5f %.5f %.5f %.5f" % (n, ps, es, pr, er)) print("BLOCK sensitivity variant short_ok short_res long_ok long_res") for nm, (sa, sb) in e6.items(): print(" %-28s %.4f %.4f %.4f %.4f" % (nm.replace(" ", "_"), sa["p_success"], sa["p_resist"], sb["p_success"], sb["p_resist"])) print("BLOCK pd C psi_S psi_R") for c in np.linspace(0.0, 32.0, 129): print(" %.4f %.6f %.6f" % (c, float(psi(c, PSI_MAX_S, MIC_S)), float(psi(c, PSI_MAX_R, MIC_R)))) print("BLOCK window MSC MIC_S MIC_R") print(" %.9f %.9f %.9f" % (msc_bis, mic_s_found, mic_r_found)) head("TOTALS") total = (n_v5 + 4000 + 3 * 2000 + N_E1 * len(doses) * len(taus) + N_E2 * len(plans) + N_E3 * 7 + N_E4 * 16 + N_E4 * 4 + N_E5 + N_E6 * 2 * len(variants) + 2) print("simulated treatment courses in this run : %d" % total) print("replicates in the largest single cell : %d" % N_E5) print("master seed : %d" % SEED) print("wall clock : %.1f s" % (time.time() - t_start)) print() print("Reminder: this is a model. It contains no patients, no bacteria and no") print("clinical data, and nothing in it is advice about taking medicine.") if __name__ == "__main__": main()