allometry_and_fih.py · Scientific Agent Skills

1#!/usr/bin/env python3

2"""Allometric scaling, maturation, and first-in-human starting dose.

3

4Two different jobs share this script because they share a failure mode: taking

5a number derived for one purpose and using it for another. An HED is not a

6starting dose. A NOAEL-derived MRSD is insufficient alone for high-risk

7immunomodulators. Fixed allometry and generic maturation need drug-specific

8justification and uncertainty assessment before extrapolation.

9

10 python3 allometry_and_fih.py --scale --cl 5 --weight-from 70 --weight-to 15

11 python3 allometry_and_fih.py --scale --cl 5 --weight-from 70 --weight-to 6 --pma-weeks 44

12 python3 allometry_and_fih.py --exponent -i species.csv

13 python3 allometry_and_fih.py --fih --noael rat=50,dog=10 --safety-factor 10 --human-weight 60

14 python3 allometry_and_fih.py --mabel --ec50 2.5 --target-occupancy 0.2 --cl 0.2

15

16``--exponent`` input is a table of ``species,weight,cl`` for cross-species

17regression.

18"""

19

20from __future__ import annotations

21

22import argparse

23import math

24from typing import Sequence

25

26import numpy as np

27from scipy.stats import t as t_dist

28

29from _common import (

30 InputError,

31 Report,

32 add_format_argument,

33 main_wrapper,

34 parse_float,

35 read_table,

36 require_columns,

37)

38

39# FDA, "Estimating the Maximum Safe Starting Dose in Initial Clinical Trials

40# for Therapeutics in Adult Healthy Volunteers" (July 2005), Table 1. Km is the

41# body-weight-to-surface-area factor; HED in mg/kg = animal mg/kg * Km_animal /

42# Km_human.

43KM_FACTORS = {

44 "human-adult": (60.0, 37.0),

45 "human-child": (20.0, 25.0),

46 "mouse": (0.020, 3.0),

47 "hamster": (0.080, 5.0),

48 "rat": (0.150, 6.0),

49 "ferret": (0.300, 7.0),

50 "guinea-pig": (0.400, 8.0),

51 "rabbit": (1.8, 12.0),

52 "dog": (10.0, 20.0),

53 "monkey": (3.0, 12.0),

54 "marmoset": (0.350, 6.0),

55 "squirrel-monkey": (0.600, 7.0),

56 "baboon": (12.0, 20.0),

57 "micro-pig": (20.0, 27.0),

58 "mini-pig": (40.0, 35.0),

59}

60

61# Anderson & Holford sigmoidal maturation on post-menstrual age. The defaults

62# are the widely used generic clearance values; a drug with a known

63# ontogeny profile should override them.

64DEFAULT_TM50_WEEKS = 54.2

65DEFAULT_MATURATION_HILL = 3.92

66

67

68def maturation_fraction(pma_weeks: float, tm50: float = DEFAULT_TM50_WEEKS, hill: float = DEFAULT_MATURATION_HILL) -> float:

69 """Fraction of adult clearance attributable to enzyme maturation."""

70 if pma_weeks <= 0:

71 raise InputError("post-menstrual age must be positive")

72 return pma_weeks**hill / (tm50**hill + pma_weeks**hill)

73

74

75def allometric(value: float, weight_from: float, weight_to: float, exponent: float) -> float:

76 return value * (weight_to / weight_from) ** exponent

77

78

79def hed_mg_per_kg(animal_dose_mg_kg: float, species: str) -> float:

80 if species not in KM_FACTORS:

81 raise InputError(f"unknown species {species!r}; known: {', '.join(sorted(KM_FACTORS))}")

82 _, km_animal = KM_FACTORS[species]

83 _, km_human = KM_FACTORS["human-adult"]

84 return animal_dose_mg_kg * km_animal / km_human

85

86

87# ---------------------------------------------------------------- exponent

88

89

90def fit_exponent(weights: np.ndarray, values: np.ndarray) -> dict[str, float]:

91 """Log-log regression of a parameter on body weight across species."""

92 if len(weights) < 3 or np.ptp(weights) <= 0:

93 raise InputError("cross-species regression needs at least 3 species and varied weights")

94 x = np.log(weights)

95 y = np.log(values)

96 n = len(x)

97 slope, intercept = np.polyfit(x, y, 1)

98 fitted = intercept + slope * x

99 ss_res = float(np.sum((y - fitted) ** 2))

100 ss_tot = float(np.sum((y - y.mean()) ** 2))

101 se_slope = math.sqrt(ss_res / (n - 2) / float(np.sum((x - x.mean()) ** 2))) if n > 2 else float("nan")

102 return {

103 "exponent": float(slope),

104 "coefficient": float(math.exp(intercept)),

105 "se_exponent": se_slope,

106 "ci95_low": float(slope - float(t_dist.ppf(0.975, n - 2)) * se_slope),

107 "ci95_high": float(slope + float(t_dist.ppf(0.975, n - 2)) * se_slope),

108 "r2": 1.0 - ss_res / ss_tot if ss_tot > 0 else float("nan"),

109 "n_species": n,

110 }

111

112

113def rule_of_exponents(exponent: float) -> str:

114 """Mahmood & Balian's rule of exponents, stated with its own caveat."""

115 if exponent < 0.55:

116 return "below 0.55: simple allometry is unreliable; human clearance is likely overpredicted"

117 if exponent <= 0.70:

118 return "0.55-0.70: simple allometry"

119 if exponent <= 1.00:

120 return "0.71-1.00: historical rule suggests lifespan scaling; do not apply without drug-specific justification"

121 return "above 1.00: historical rule suggests brain-weight scaling; do not apply without drug-specific justification"

122

123

124# --------------------------------------------------------------------- CLI

125

126

127def build_parser() -> argparse.ArgumentParser:

128 parser = argparse.ArgumentParser(

129 description="Allometric scaling, paediatric maturation, and first-in-human dose estimation.",

130 formatter_class=argparse.RawDescriptionHelpFormatter,

131 epilog=__doc__,

132 )

133 mode = parser.add_argument_group("mode (choose one)")

134 mode.add_argument("--scale", action="store_true", help="scale parameters between two body weights")

135 mode.add_argument("--exponent", action="store_true", help="estimate the allometric exponent from species data")

136 mode.add_argument("--fih", action="store_true", help="NOAEL-based human equivalent dose and MRSD")

137 mode.add_argument("--mabel", action="store_true", help="illustrative equilibrium occupancy/dose-rate calculation; not a complete MABEL")

138

139 parser.add_argument("-i", "--input", help="species table for --exponent (species,weight,cl)")

140 parser.add_argument("--cl", type=float, help="clearance at the reference weight")

141 parser.add_argument("--volume", type=float, help="volume at the reference weight")

142 parser.add_argument("--weight-from", type=float, default=70.0, help="reference weight, kg (default: 70)")

143 parser.add_argument("--weight-to", type=float, help="target weight, kg")

144 parser.add_argument("--cl-exponent", type=float, default=0.75)

145 parser.add_argument("--v-exponent", type=float, default=1.0)

146 parser.add_argument("--pma-weeks", type=float, help="post-menstrual age; adds a maturation factor")

147 parser.add_argument("--tm50", type=float, default=DEFAULT_TM50_WEEKS)

148 parser.add_argument("--maturation-hill", type=float, default=DEFAULT_MATURATION_HILL)

149

150 parser.add_argument("--noael", help="comma-separated species=mg/kg, e.g. rat=50,dog=10")

151 parser.add_argument("--safety-factor", type=float, default=10.0)

152 parser.add_argument("--human-weight", type=float, default=60.0, help="kg, for converting MRSD to a total dose")

153

154 parser.add_argument("--ec50", type=float, help="in vitro potency (same units as the target concentration)")

155 parser.add_argument("--target-occupancy", type=float, default=0.10, help="fraction of target engagement to aim for")

156 parser.add_argument("--weight", type=float, default=70.0)

157 parser.add_argument("--mabel-safety-factor", type=float, default=1.0)

158 add_format_argument(parser)

159 return parser

160

161

162def run(argv: Sequence[str] | None = None) -> int:

163 args = build_parser().parse_args(argv)

164 from _common import validate_numeric_args

165 validate_numeric_args(args)

166 chosen = [m for m in ("scale", "exponent", "fih", "mabel") if getattr(args, m)]

167 if len(chosen) != 1:

168 raise InputError("choose exactly one of --scale, --exponent, --fih, --mabel")

169 report = Report()

170 for name in ("weight_from", "weight_to", "human_weight", "weight", "safety_factor", "mabel_safety_factor", "tm50", "maturation_hill", "cl", "volume", "ec50", "pma_weeks"):

171 value = getattr(args, name)

172 if value is not None and value <= 0:

173 raise InputError(f"--{name.replace('_', '-')} must be positive")

174

175 if args.scale:

176 if args.weight_to is None:

177 raise InputError("--scale needs --weight-to")

178 if args.cl is None and args.volume is None:

179 raise InputError("--scale needs --cl and/or --volume")

180 rows = []

181 maturation = maturation_fraction(args.pma_weeks, args.tm50, args.maturation_hill) if args.pma_weeks else 1.0

182 if args.cl is not None:

183 size_only = allometric(args.cl, args.weight_from, args.weight_to, args.cl_exponent)

184 rows.append(

185 {

186 "parameter": "CL",

187 "reference": args.cl,

188 "exponent": args.cl_exponent,

189 "size_scaled": size_only,

190 "maturation_factor": maturation,

191 "final": size_only * maturation,

192 }

193 )

194 if args.volume is not None:

195 size_only = allometric(args.volume, args.weight_from, args.weight_to, args.v_exponent)

196 rows.append(

197 {

198 "parameter": "V",

199 "reference": args.volume,

200 "exponent": args.v_exponent,

201 "size_scaled": size_only,

202 "maturation_factor": 1.0,

203 "final": size_only,

204 }

205 )

206 report.table(f"scaled from {args.weight_from} kg to {args.weight_to} kg", rows)

207 report.note("This helper applies maturation to clearance only; developmental volume changes need a separate justified model.")

208 if args.pma_weeks:

209 report.scalar("post_menstrual_age_weeks", args.pma_weeks)

210 report.scalar("maturation_fraction_of_adult_cl", maturation)

211 report.note(

212 f"maturation uses TM50 = {args.tm50} weeks PMA and Hill = {args.maturation_hill} "

213 "(illustrative). Justify drug/pathway-specific ontogeny and uncertainty before extrapolation."

214 )

215 else:

216 if args.weight_to < 20:

217 report.finding(

218 "scaling below 20 kg using size alone; assess age, drug-specific maturation and organ function. A generic PMA multiplier does not validate the prediction."

219 )

220 report.note("allometry is a covariate model for size, not evidence of a mechanism")

221

222 elif args.exponent:

223 if not args.input:

224 raise InputError("--exponent needs -i with columns species,weight,cl")

225 rows = read_table(args.input)

226 require_columns(rows, ["weight", "cl"], str(args.input))

227 weights = np.asarray([parse_float(r["weight"], "weight") for r in rows], dtype=float)

228 values = np.asarray([parse_float(r["cl"], "cl") for r in rows], dtype=float)

229 if np.any(weights <= 0) or np.any(values <= 0):

230 raise InputError("weights and clearances must be positive")

231 fit = fit_exponent(weights, values)

232 for key, value in fit.items():

233 report.scalar(key, value)

234 report.scalar("rule_of_exponents", rule_of_exponents(fit["exponent"]))

235 report.table(

236 "observed vs fitted",

237 [

238 {

239 "species": r.get("species", f"#{i + 1}"),

240 "weight_kg": w,

241 "cl_observed": v,

242 "cl_fitted": fit["coefficient"] * w ** fit["exponent"],

243 "fold_error": (fit["coefficient"] * w ** fit["exponent"]) / v,

244 }

245 for i, (r, w, v) in enumerate(zip(rows, weights, values))

246 ],

247 )

248 if fit["n_species"] < 4:

249 report.finding(

250 f"exponent estimated from {fit['n_species']} species; the confidence interval "

251 f"({fit['ci95_low']:.2f} to {fit['ci95_high']:.2f}) is too wide to distinguish 0.75 from "

252 "most alternatives"

253 )

254 if not (fit["ci95_low"] <= 0.75 <= fit["ci95_high"]):

255 report.note("the 95% interval excludes 0.75; a fixed-exponent model would be misspecified here")

256 report.note(

257 "Cross-species allometry depends on species relevance, binding, pathways and input uncertainty. It is one exploratory input, not a validated starting-dose prediction."

258 )

259

260 elif args.fih:

261 if not args.noael:

262 raise InputError("--fih needs --noael, e.g. --noael rat=50,dog=10")

263 entries = []

264 for chunk in args.noael.split(","):

265 if "=" not in chunk:

266 raise InputError(f"--noael entries look like species=mg/kg; got {chunk!r}")

267 species, dose_text = chunk.split("=", 1)

268 species = species.strip().lower()

269 dose = parse_float(dose_text, f"NOAEL for {species}")

270 if dose is None or dose <= 0:

271 raise InputError(f"NOAEL for {species} must be positive")

272 hed = hed_mg_per_kg(dose, species)

273 weight, km = KM_FACTORS[species]

274 entries.append(

275 {

276 "species": species,

277 "noael_mg_kg": dose,

278 "reference_weight_kg": weight,

279 "km": km,

280 "divide_by": KM_FACTORS["human-adult"][1] / km,

281 "hed_mg_kg": hed,

282 }

283 )

284 report.table("human equivalent dose by species", entries)

285 most_sensitive = min(entries, key=lambda e: e["hed_mg_kg"])

286 mrsd = most_sensitive["hed_mg_kg"] / args.safety_factor

287 report.scalar("lowest_hed_species", most_sensitive["species"])

288 report.scalar("lowest_hed_mg_kg", most_sensitive["hed_mg_kg"])

289 report.scalar("safety_factor", args.safety_factor)

290 report.scalar("mrsd_mg_kg", mrsd)

291 report.scalar("mrsd_total_mg", mrsd * args.human_weight)

292 report.note(

293 "HED conversion uses body-surface-area scaling from FDA's 2005 maximum-safe-starting-dose "

294 "guidance, Table 1. It applies to small molecules; for biologics whose clearance is not "

295 "surface-area-related, mg/kg or exposure matching is generally more appropriate."

296 )

297 report.note(

298 "The calculator selects the lowest numerical HED. Establish pharmacological/toxicological species relevance before choosing a starting-dose basis."

299 )

300 if args.safety_factor < 10:

301 report.finding(

302 f"safety factor of {args.safety_factor} is below the default 10; a reduction has to be "

303 "justified, and an increase is expected for steep dose-response, irreversible toxicity, "

304 "novel targets, or nonlinear PK"

305 )

306 report.finding(

307 "MRSD from a NOAEL is not appropriate on its own for agonist immunomodulators or other "

308 "agents with a plausible risk of severe on-target toxicity. Integrate NOAEL, pharmacologically active dose, MABEL, exposure and uncertainty; --mabel is only a simplified calculation"

309 )

310

311 else: # --mabel

312 if args.ec50 is None:

313 raise InputError("--mabel needs --ec50")

314 if not 0 < args.target_occupancy < 1:

315 raise InputError("--target-occupancy must be between 0 and 1")

316 if args.cl is None:

317 raise InputError("--mabel needs --cl to calculate a model-implied dose rate")

318 target_conc = args.ec50 * args.target_occupancy / (1.0 - args.target_occupancy)

319 # Steady-state-equivalent dose to reach the target average concentration.

320 dose_rate = target_conc * args.cl

321 report.scalar("ec50", args.ec50)

322 report.scalar("target_fractional_effect_or_occupancy", args.target_occupancy)

323 report.scalar("target_concentration", target_conc)

324 report.scalar("clearance", args.cl)

325 report.scalar("dose_rate_for_target_concentration", dose_rate)

326 report.scalar("dose_rate_with_safety_factor", dose_rate / args.mabel_safety_factor)

327 report.note(

328 "This simplified calculation inverts a fractional Emax/equilibrium occupancy curve. A full MABEL assessment integrates human biology, exposure, target engagement and uncertainty."

329 )

330 report.note(

331 "The EC50 formula describes a fractional Emax response. For receptor occupancy, supply a justified Kd instead: functional EC50 need not equal Kd. Concentration times CL is amount/time (e.g. mg/L times L/h gives mg/h), not an initial dose. TMDD and time-varying target require a mechanistic model."

332 )

333 report.finding(

334 "MABEL requires human-cell in vitro potency, target expression in the relevant tissue, and "

335 "a defined concentration-occupancy relationship. Confirm all three are drug-specific before "

336 "using this number."

337 )

338

339 return report.emit(args.format)

340

341

342if __name__ == "__main__":

343 raise SystemExit(main_wrapper(run))