""" PRĀṆA–PD v2.0 — Python reference implementation (numpy only). dx/dt = [A_L + rho(t) A_C - D(t)] x + h + B u(t) + G E(t) + w(t) Mirrors src/model.ts (the web simulator) so the two can be cross-validated: python python/pranapd.py --validate All coefficients are ILLUSTRATIVE PLACEHOLDERS to be estimated from data. Research use only: not a medical device, not a dosing or prescribing tool. """ import argparse, json, math import numpy as np DOMAINS = ["Prana", "Apana", "Samana", "Vyana", "Udana"] PROFILES = {"A": [.20, .15, .20, .75, .25], # predominantly motor "B": [.15, .75, .70, .30, .15], # autonomic / GI "C": [.35, .10, .25, .25, .80]} # speech / swallowing LAMBDA = np.array([.16, .14, .14, .18, .16]) THETA = np.array([.25, .25, .25, .25, .30]) R_DIAG = np.array([.02, .05, .03, .03, .04, .02, .04]) G_DIAG = -0.12 EVENT_SHAPE = np.array([1.0, .3, .5, .6, .8]) BMAT = np.array([[0, .015, 0, .06, .01, .02, .01], [0, .02, .06, 0, 0, .01, .01], [.01, .01, .06, 0, 0, .01, .01], [.04, .07, .01, .01, 0, .015, .005], [.01, .005, .005, .02, .07, .005, 0]]) PAIRS = {(0, 4): .10, (1, 2): .10, (2, 3): .07, (3, 4): .07, (0, 3): .07} W = np.zeros((5, 5)) for (i, j), v in PAIRS.items(): W[i, j] = W[j, i] = v def laplacian(C=None): """A_C: weighted graph Laplacian (zero row sums) of the masked coupling weights.""" C = np.ones((5, 5)) if C is None else np.asarray(C) M = W * np.minimum(C, C.T) np.fill_diagonal(M, 0) return M - np.diag(M.sum(1)) def d_from_s(s, lam=LAMBDA): s = np.clip(s, 0, .95) return lam * s / (1 - s) def simulate(profile="A", arm="usual", weeks=52, rho0=.8, rho_decay=.0005, progression=None, goal=.5, kappa=3., review=4, q=None, lam_safe=2., adjunct=1., event=None, usual=.5, sigma_w=0., sigma_m=0., seed=7, dt=.1, coupling=None, s0=None): s0 = np.array(PROFILES[profile] if s0 is None else s0, float) prog = np.array([.001, .001, .001, .0015, .001] if progression is None else progression) q = np.ones(5) if q is None else np.asarray(q) AC = laplacian(coupling) Bm = BMAT.copy() Bm[:, 6] *= adjunct if arm == "adaptive_ayur" else 0. x0 = 1 - np.clip(s0, 0, .95) xstar = x0 + goal * (1 - x0) K = np.zeros((7, 5)) for j in range(1, 7): if j == 6 and arm != "adaptive_ayur": continue K[j] = kappa * BMAT[:, j] / .06 rng_m, rng_w = np.random.default_rng(seed), np.random.default_rng(seed + 1) x = x0.copy(); u = np.array([usual, 0, 0, 0, 0, 0, 0.]) steps = int(round(weeks / dt)); rsteps = int(round(review / dt)) traj, Jt, Je, Js = [], 0., 0., 0. for k in range(steps + 1): t = k * dt rho = float(np.clip(rho0 - rho_decay * t, 0, 1)) if arm != "usual" and k % rsteps == 0: xhat = np.clip(x + sigma_m * rng_m.standard_normal(5), 0, 1) e = xstar - xhat u = np.clip(K @ e, 0, 1); u[0] = usual if k % 10 == 0: e = xstar - x S = np.sum((np.maximum(0, THETA - x) / THETA) ** 2) traj.append(dict(t=int(round(t)), x=x.copy(), u=u.copy(), rho=rho)) if round(t) >= 1: Jt += float(q @ e ** 2); Je += float(R_DIAG @ u ** 2); Js += lam_safe * S if k == steps: break s = np.clip(s0 + prog * t, 0, .95) E = (EVENT_SHAPE * event["mag"] if event and event["week"] <= t < event["week"] + event.get("dur", 3) else 0) dx = -LAMBDA * x - d_from_s(s) * x + rho * AC @ x + LAMBDA + Bm @ u + G_DIAG * E x = np.clip(x + dt * dx + math.sqrt(dt) * sigma_w * rng_w.standard_normal(5), 0, 1) return dict(traj=traj, xstar=xstar, J=Jt + Je + Js, Jtrack=Jt, Jeffort=Je, Jsafe=Js, final=traj[-1]["x"], mean=float(traj[-1]["x"].mean())) def steady_state(s, rho, u, lam=LAMBDA, coupling=None): AC = laplacian(coupling) M = np.diag(lam + d_from_s(np.asarray(s))) - rho * AC return np.clip(np.linalg.solve(M, lam + BMAT @ np.asarray(u)), 0, 1) def eigen_report(s, rho, coupling=None): M = -np.diag(LAMBDA + d_from_s(np.asarray(s))) + rho * laplacian(coupling) ev = np.linalg.eigvalsh((M + M.T) / 2) return dict(max_real=float(ev.max()), stable=bool(ev.max() < 0)) # ------------------------------------------------------------------ identification def synth_cohort(n=60, weeks=52, sigma_w=.02, sigma_m=.015, seed=1, rho=.8): """Simulate a cohort with KNOWN coupling and weekly noisy measurements. Random, RECORDED intercurrent events (e.g. infection, medication change) provide the excitation that identification needs. Returns (list of X[t, 5], list of event indicators).""" rng = np.random.default_rng(seed) out, evs = [], [] for p in range(n): s0 = np.clip(rng.beta(2, 5, 5) * 1.1, .03, .85) w0, mag = rng.uniform(4, 40), rng.uniform(.5, 2) sim = simulate(s0=s0, arm="usual", weeks=weeks, rho0=rho, rho_decay=0, event=dict(week=w0, mag=mag, dur=3), sigma_w=sigma_w, seed=int(rng.integers(1e6)), usual=0.) X = np.array([r["x"] for r in sim["traj"]]) out.append(X + sigma_m * rng.standard_normal(X.shape)) tm = np.arange(len(X) - 1) + .5 # interval midpoints evs.append(mag * ((tm >= w0) & (tm < w0 + 3)).astype(float)) return out, evs def estimate_coupling(cohort, events=None, dt=1., pairs=tuple(PAIRS), iv=True): """Identify pairwise coupling rho*w_ij from repeated noisy measurements. Per domain i (weekly observations, trapezoid rule on the interval [t, t+1]): x_i(t+1) - x_i(t) = dt * [ c_i - a_{i,p} xbar_i + sum_j g_ij (xbar_j - xbar_i) ] where xbar = (x(t) + x(t+1)) / 2. * a_{i,p} (= lambda_i + d_i) is PATIENT-SPECIFIC; forcing a shared value biases the coupling toward zero. g_ij is shared across patients. * Measurement noise enters the regressors AND the response, biasing OLS upward (errors-in-variables). With iv=True the regressors are instrumented by the observations two time-points earlier, x(t-1), whose noise is independent of the noise in x(t), x(t+1) (two-stage least squares). * Only structurally hypothesised pairs enter (structural zeros of A_C). * RECORDED disturbances E(t) enter as a regressor (shared sensitivity g_i). Omitting a common shock makes coupling absorb it (found during validation: the biggest bias). """ n = len(cohort) est = {} for i in range(5): nbrs = [j for j in range(5) if (i, j) in pairs or (j, i) in pairs] F, Z, y = [], [], [] for p, X in enumerate(cohort): T = len(X) - 2 # intervals t = 1 .. len-2 (need x(t-1)) xb = (X[2:] + X[1:-1]) / 2 # midpoint of [t, t+1] xl = X[:-2] # x(t-1): instrument for x(t) resp = (X[2:, i] - X[1:-1, i]) / dt own = np.zeros((T, n)); own[:, p] = -xb[:, i] own_z = np.zeros((T, n)); own_z[:, p] = -xl[:, i] cpl = np.column_stack([xb[:, j] - xb[:, i] for j in nbrs]) cpl_z = np.column_stack([xl[:, j] - xl[:, i] for j in nbrs]) ecol = np.zeros((T, 1)) if events is None else events[p][1:T + 1].reshape(-1, 1) # interval t (t=1..T) F.append(np.column_stack([np.ones(T), own, cpl, ecol])) Z.append(np.column_stack([np.ones(T), own_z, cpl_z, ecol])) y.append(resp) F, Z, y = np.vstack(F), np.vstack(Z), np.concatenate(y) if iv: Pz = Z @ np.linalg.pinv(Z.T @ Z) @ Z.T Fh = Pz @ F beta = np.linalg.lstsq(Fh.T @ F, Fh.T @ y, rcond=None)[0] else: beta = np.linalg.lstsq(F, y, rcond=None)[0] for k, j in enumerate(nbrs): est.setdefault((min(i, j), max(i, j)), []).append(beta[1 + n + k]) return {k: float(np.mean(v)) for k, v in est.items()} def bootstrap(cohort, events=None, B=40, seed=0): rng = np.random.default_rng(seed); draws = [] for _ in range(B): idx = rng.integers(0, len(cohort), len(cohort)) draws.append(estimate_coupling([cohort[i] for i in idx], None if events is None else [events[i] for i in idx])) return {k: (float(np.percentile([d[k] for d in draws], 2.5)), float(np.percentile([d[k] for d in draws], 97.5))) for k in draws[0]} def main(): ap = argparse.ArgumentParser() ap.add_argument("--validate", action="store_true", help="print deterministic summaries for cross-check with the TypeScript model") ap.add_argument("--identify", action="store_true", help="run the coupling-recovery experiment") a = ap.parse_args() if a.validate: rows = {} for pr in "ABC": for arm in ("usual", "adaptive", "adaptive_ayur"): r = simulate(pr, arm) rows[f"{pr}/{arm}"] = dict(mean=round(r["mean"], 4), J=round(r["J"], 4), final=[round(float(v), 2) for v in r["final"]]) print(json.dumps(rows, indent=1)) if a.identify: coh, ev = synth_cohort() est = estimate_coupling(coh, ev); ci = bootstrap(coh, ev) print("pair true(rho*w) estimate 95% bootstrap CI") for (i, j), w in PAIRS.items(): e = est[(i, j)]; lo, hi = ci[(i, j)] print(f"{DOMAINS[i]:7s}-{DOMAINS[j]:7s} {0.8*w:8.3f} {e:8.3f} [{lo:6.3f}, {hi:6.3f}]") if __name__ == "__main__": main()