APOD: reasoning-guided agentic population ordinary differential equation discovery for pharmacological digital twins

1 Introduction

• Automated ODE methods either analyze one subject at a time or optimize a numeric fit-and-parsimony objective without accounting for biological plausibility.Biomedical datasets often contain sparse, noisy, partial observations of mechanisms shared across subjects.Source: §1 · p. 2
• Existing LLM agents lack population-level inference, allow only predefined model dimensions, rely on a single fit metric, and have rarely been tested on sparse heterogeneous cohorts.Automated NLME discovery is also limited to predefined pharmacokinetic building blocks, while neural-ODE approaches require dense subject sampling and large cohorts.Source: §1 · p. 2
• A population digital twin uses a shared ODE structure with subject-specific parameters drawn from a population distribution, supporting personalized prediction and counterfactual simulation.Sparse measurements can personalize the population fit to forecast an endpoint such as treatment-toxicity probability or simulate an alternative dosing regimen.Source: §1 · p. 2
• APOD is introduced to combine open structural search with population inference, and is evaluated on synthetic benchmarks, real experimental cohorts, and radioligand-therapy platelet dynamics.The platelet model is assessed for held-out-cohort forecasts.Source: §1 · p. 2

2.1 APOD: Agentic Population ODE Discovery

• APOD proposes both a shared ODE and between-subject variability model, then iteratively evaluates fit, identifiability, and plausibility to select a model.The selected model is returned with its diagnostics and reasoning trace; the prompts encode no domain-specific knowledge.Source: §2.1 · pp. 2–3

2.2 APOD recovers ground-truth ODEs by reasoning over fit diagnostics

• 96% overall structural recovery; APOD recovered PK structures in 94% of runs versus mlxModelFinder’s 98%, and TD structures in 100% versus 80%.APOD’s median run was 24 times faster on PK datasets and 4.5 times faster on TD scenarios.Source: §2.2 · p. 3
• Without reasoning, structural recovery fell from 96% to 88%, exact random-effect recovery from 70% to 59%, and median runtime rose from 4.1 to 10.8 min.With reasoning, later iterations corrected the first proposal in 44 of the 49 runs where it was wrong.Source: §2.2 · p. 5

2.3 APOD suggests relevant structures on real clinical and experimental cohorts

• From a drug-agnostic prompt, APOD selected a three-compartment IV docetaxel model with identifiable parameters and prediction intervals covering observed percentiles in training and held-out patients.It also reached plausible, identifiable disposition models on oral ciclosporine and [177Lu]Lu-DOTATATE renal-handling cohorts.Source: §2.3 · p. 7
• On real tumor cohorts, APOD selected hyper-Gompertz untreated growth and logistic growth with a threshold drug effect through an effect compartment, with good held-out predictions.For untreated growth, it removed variability from two parameters after diagnosing high shrinkage; the treated-tumor model predicted individual and population behavior across dose levels and schedules.Source: §2.3 · pp. 7–8

2.4 APOD provides a PDT that predicts thrombocytopenia under radioligand therapy

• For radioligand-therapy platelet dynamics, APOD introduced a dose-driven drug proxy and recovered a seven-ODE system with nonlinear feedback and identifiable parameters.The model described the cohort’s platelet-toxicity time course without signs of misspecification in individual fits or visual predictive checks.Source: §2.4 · p. 9
• In 19 later-period patients, first-cycle individualisation from a median of three measurements produced 90% predictive intervals containing the observed nadir in all patients.Predicted nadirs had concordance correlation 0.91; the population model without individualisation explained almost none of the between-patient nadir variation (R2 0.06).Source: §2.4 · p. 9
• For a patient predicted to cross the grade-1 threshold, halving injected activity from cycle two lowered predicted grade-1 nadir probability from 0.94 to 0.59, while omitting a cycle left it at 0.93.Across 19 patients, predicted nadirs rose by 26 G/L with halved activity, 7 G/L with an omitted cycle, and 2 G/L with a four-week delay.Source: §2.4 · p. 9

3 Discussion

• APOD matched or exceeded structural recovery 2.5- to 40-fold faster than a library-based tool, fitting about six times fewer candidate models.The comparison was with an automated model-selection tool from the pharmacometrics community.Source: §3 · p. 10
• On synthetic benchmarks, APOD recovered the exact data-generating ODE in most runs within minutes under sparse sampling, few subjects, and unobserved latent compartments.Reasoning across fit, identifiability, and plausibility outperformed ranking candidates on BICc alone: recovery was 96% versus 88%, and single-metric runs took two to three times longer.Source: §3 · pp. 10–11
• The APOD platelet model was discovered de novo from a blinded one-sentence prompt; first-cycle forecasts for three test patients had mean absolute relative errors of 4%, 8%, and 22%, respectively.For patient 1, cycle-1 data were used for individualisation before simulation under the standard 7,400 MBq schedule, one omitted cycle, or 3,700 MBq from cycle 2.Source: Figure 3 · p. 10
• Structural recovery was more consistent than recovery of between-subject variability, especially when simulated variance was small relative to residual error, limiting confidence in individual-prediction widths.The authors identify variability recovery as a clear area for improvement.Source: §3 · p. 11
• From one measured platelet count, APOD discovered a seven-ODE Friberg-type model de novo, and first-cycle forecasts contained observed platelet nadirs in all 19 patients of a temporally independent cohort.The test cohort comprised 19 patients; the figure shows first-cycle forecast errors of 4%, 8%, and 22% for three test patients.Source: Figure 3 · p. 10, §3 · p. 11
• APOD does not model covariates, and severe dose-limiting-event prediction remains unestablished because clinical validation was temporal rather than multi-centre and observed toxicity was moderate.Discovery in domains without an existing mechanistic description also remains untested.Source: §3 · p. 11

4 Methods

• APOD can run with cloud-hosted language models through an API or local open-weight models through Ollama without changing the discovery loop.All experiments in the paper used Claude Opus 4.5, and the seven-backbone benchmark changed only the language model while keeping the loop, prompts, stopping rules, and estimation settings identical.Source: §4.1 · p. 12
• The language model received dataset structure and aggregated initialization metrics, never individual-level records, while individual data remained local to Monolix.The summaries included cohort baselines, overall means, and late-phase values with their mean times, which anchored plausible initial parameter values.Source: §4.1 · p. 12
• Structure-agnostic prompts named the observed quantity and task but supplied no mechanism or model family.The platelet prompt asked for a platelet-count model under radioligand therapy; PK and tumor-dynamics prompts were similarly task-based, with wording adapted to the dataset.Source: §4.1 · p. 12
• The discovery loop cycles from a prompt and aggregated metrics through candidate generation, local fitting, and evaluation to structural refinement or stopping.A worked PK example describes the agent proposing initial models, writing valid Mlxtran code, and receiving a structured judgment on fit, identifiability, and plausibility.Source: §4.1 · p. 12

4.2 Model representation, estimation and iterative refinement

• Candidate models were explicit Mlxtran ODE systems whose equations and structures were agent-generated rather than restricted to a fixed template library.The output variable and residual-error model were inputs, specified from the data-generating model for synthetic benchmarks or from published analyses or analyst inspection for real datasets.Source: §4.2 · p. 12
• Population parameters were estimated by SAEM, with log-normal inter-individual variability on structural parameters that the LLM could remove based on diagnostics.Runs used one Markov chain with automatic stopping disabled; the standard benchmark used 1,000 exploratory and 500 smoothing iterations, with a more complex PK-TD case using 1,500 exploratory iterations.Source: §4.2 · p. 13
• After each fit, the agent received fit, identifiability, residual, and predictive-check diagnostics and reasoned over them to choose the next structural move or stop.Each run allowed up to ten structural iterations; the agent weighed model quality, convergence, and parameter plausibility against further exploration.Source: §4.2 · p. 13

4.3 Ablation of the reasoning layer

• The reasoning ablation reduced feedback to objective function, information criteria, and parameter count, removed reasoning and stopping instructions, and selected the lowest-BICc candidate.Five runs were performed per dataset; the agent could still stop by proposing no further candidate.Source: §4.3 · p. 13

4.4 Automated model selection in population modeling with mlxModelFinder

• mlxModelFinder searched standard PK and curated tumor-growth libraries using the same SAEM settings and the same machine to compare structural search and runtime.PK searches used clearance parameterisation and five estimation seeds under a 12-hour budget; one of 50 searches exceeded that limit.Source: §4.4 · p. 13
• For PK scenarios, mlxModelFinder used ant-colony optimisation to search inter-individual-variability structure while holding the residual-error model fixed and scoring candidates with a gradient-boosting surrogate.The search used population size 20 and up to ten iterations.Source: §4.4 · p. 14

4.5 Synthetic benchmarks

• Ten synthetic pharmacokinetic scenarios span one- to three-compartment intravenous and oral models, while five tumor-growth structures generate untreated baseline datasets.Datasets were designed for practical parameter identifiability, checked using SAEM initialized at ground-truth values; APOD was run ten times per dataset.Source: §4.5 · p. 14

4.6.1 Real-world model discovery: PK, TD and drug–response coupling

• Real-world datasets cover docetaxel, ciclosporine, radioligand renal handling, two untreated tumor models, and paclitaxel exposure coupled to tumor volumes across routes and schedules.For the paclitaxel task, previously validated population PK parameters generated individual exposures; discovery used those exposures and measured tumor volumes.Source: §4.6.1 · pp. 14–15

4.6.2 Population digital twin construction in clinical radioligand therapy

• The COMPORIV cohort provides 519 platelet measurements from 44 radioligand-therapy patients, with 19 later-treated patients reserved for temporal validation using first-cycle data.The test cohort supports patient-specific platelet-trajectory and nadir predictions without covariates.Source: §4.6.2 · p. 15

4.7 Evaluation of APOD’s discovery and descriptive abilities

• Synthetic runs count as structural recoveries when the selected model reproduces the known data-generating structure.The criterion and per-run reporting are specified under Statistics.Source: §4.7 · p. 15
• Real-cohort descriptive adequacy is judged from individual profiles and observed-versus-predicted plots, plus population visual predictive checks.A model is considered adequate when observed percentiles fall within the prediction intervals across the sampling window.Source: §4.7 · p. 15
• Practical identifiability on real cohorts is assessed using population-parameter relative standard errors and random-effect eta-shrinkage.Source: §4.7 · p. 15

4.8 Held-out evaluation on real datasets

• Real PK and tumor-dynamics cohorts with at least 20 individuals were randomly split 70/30, keeping each individual’s observations on one side.The split used a fixed seed, and discovery retained the full-cohort configuration.Source: §4.8 · p. 15
• For held-out individuals, typical values, between-subject variances, and residual error were frozen; only their individual random effects were estimated.Source: §4.8 · p. 15
• Ciclosporine and [177Lu]Lu-DOTATATE renal-handling cohorts were reported on the full cohort because their 14 and 10 individuals would leave only three to four held-out patients.That held-out size was considered too small for interpretable visual predictive checks or prediction-error estimates.Source: §4.8 · pp. 15–16

4.9 Individual forecasting and predictive validation

• First-cycle data individualize the fixed COMPORIV model for test patients; posterior draws predict nadirs over cycles 2–4 and simulate delayed, omitted, or halved later cycles.The counterfactuals retain the fitted structure, population parameters, and individual draws; these simulations were not compared with observations.Source: §4.9 · p. 16

4.10 Statistics

• Nadir prediction is evaluated with R2, Lin’s concordance, Spearman correlation, RMSE, MAE, and 90% interval coverage; grade-1 flagging uses sensitivity and specificity.Grade 1 is defined here as a platelet nadir below 150 G/L; population-only signal uses R2 without individualization.Source: §4.10 · p. 16

4.11 Computational environment

• Analyses ran on a single Apple MacBook Pro with an M4 Pro chip, using Python for the agent and forecasting and Monolix for population estimation and simulation.The specified environment was macOS 26.5 with 14 CPU cores and 48 GB unified memory.Source: §4.11 · p. 16

4.12 Ethics and data governance

• The COMPORIV training cohort was a sub-cohort of the retrospective, single-centre PIXIRIV study of patients with neuroendocrine tumors treated at Institut Paoli-Calmettes.The independent test cohort comprised 19 patients treated at the same centre over a later, non-overlapping period.Source: §4.12 · pp. 16–17
• Individual patient data were pseudonymised, processed in a private institutional repository under GDPR, and never exposed to the language model.Source: §4.12 · p. 17

4.13 Use of large language models

• Claude assisted with manuscript translation and language editing and portions of code development; the authors reviewed all outputs.Source: §4.13 · p. 17

4.14 Data availability

• Synthetic benchmark datasets and discovered Mlxtran model files are available from the paper’s GitLab repository.Source: §4.14 · p. 17
• Preclinical tumor-growth data are available through Zenodo, while individual clinical cohort data are available on reasonable request subject to data-sharing agreement and institutional approval.The clinical cohorts cannot be shared publicly under GDPR.Source: §4.14 · p. 17

4.15 Code availability

• The APOD code and the analysis and figure-generation code are available in separate GitLab repositories.Source: §4.15 · p. 17

• The workflow figure traces prompt and aggregated metrics through model suggestion, code writing, fitting, diagnostic evaluation, iterative refinement, and stopping.Source: APOD supplementary figures and tables · p. 21
• The synthetic-fit figure covers individual fits and visual predictive checks for ten PK and five tumor-dynamics archetypes.Its legend identifies one panel for individual fits and another for the corresponding predictive checks.Source: APOD supplementary figures and tables · p. 23
• The ciclosporine supplement reports an APOD-discovered oral PK model and parameter estimates for 14 patients and 154 observations; the cohort was not split because fewer than 30 patients made a four-patient holdout too small for interpretable checks.Source: APOD supplementary figures and tables · p. 24
• The platelet-cohort supplement describes a Friberg-type model fitted to 44 patients and 519 observations, with a proliferating pool, five maturation compartments, feedback, and linear drug inhibition.The accompanying panels show parameter estimates, individual fits, and a visual predictive check.Source: APOD supplementary figures and tables · p. 27
• For 19 independent test-cohort patients, forecasts condition on the first treatment cycle and show observed counts, posterior-median predictions, and 90% predictive intervals.Source: APOD supplementary figures and tables · p. 28
• The synthetic PK design uses ten archetypes with 60 individuals each, while the TD design uses five growth models with 30 individuals and measurements every two days over 30 days.Source: APOD supplementary figures and tables · pp. 30, 32
• The real-cohort parameter table covers docetaxel PK, xenograft tumor dynamics, coupled PK–TD response, and platelet dynamics, with estimates and uncertainty summaries.Source: APOD supplementary figures and tables · p. 34
• The platelet-cohort characteristics table describes treatment activity and cycle spacing, platelet nadir and toxicity grades; no patient in either cohort reached grade 3.Source: APOD supplementary figures and tables · p. 35