LLM-Guided ODE Discovery and Parameter Inference from Small-Cohort Aggregate Data

# LLM-Guided ODE Discovery and Parameter Inference from Small-Cohort Aggregate Data
Source: [https://arxiv.org/html/2607.00733](https://arxiv.org/html/2607.00733)
Hanning Yang1∗,Meropi Karakioulaki2,Lennart Purucker3,Tim Litwin1,Cristina Has2,Moritz Hess11Institute of Medical Biometry and Statistics, Faculty of Medicine and Medical Center, University of Freiburg, Germany 2Department of Dermatology, Medical Faculty and Medical Center, University of Freiburg, Germany 3Prior Labs, University of Freiburg, Germany hanning\.yang@uniklinik\-freiburg\.de

###### Abstract

Mechanistic modeling via ordinary differential equations \(ODEs\) provides interpretable descriptions of complex dynamics and enables inference of underlying mechanisms, which is particularly valuable in clinical settings\. However, in rare diseases, both the structure and parameters of the model are typically unknown, while individual\-level data is scarce, noisy, heterogeneous, and subject to privacy constraints\. In such settings, population\-level summary statistics provide a practical privacy\-preserving data representation, while capturing heterogeneity further requires modeling parameters as distributions rather than fixed values\. Yet no existing method jointly discovers ODE structure and refines parameter distributions solely from summary statistics\. We present AgentODE, an end\-to\-end framework that addresses this gap\. An LLM proposes candidate ODE structures, while a tool\-augmented inference agent iteratively refines parameter distributions through a diagnosis–update loop, operating on population\-level summary statistics alone\. We evaluate AgentODE on three benchmark problems across different fields and two clinical datasets, including the rare disease recessive dystrophic epidermolysis bullosa \(RDEB\), with only 231 observations across 46 patients\. AgentODE recovers functionally consistent ODE structures across all settings, and experiments on RDEB demonstrates that in sparse and noisy data settings reasoning from summary statistics promotes mechanistically principled structure discovery, whereas baselines with individual\-level data access recover implausible structures despite better predictive performance\. AgentODE opens new possibilities for mechanistic modeling of rare diseases directly from population\-level summary statistics, where data scarcity and privacy constraints have traditionally limited such analyses111Code available at[https://github\.com/HanningYang/AgentODE](https://github.com/HanningYang/AgentODE)\.\.

## 1Introduction

Mechanistic models based on ordinary differential equations \(ODEs\) provide a principled framework for describing the dynamics of complex systems across domains such as biology, medicine, and engineering\(Raueet al\.,[2013](https://arxiv.org/html/2607.00733#bib.bib14)\)\. By explicitly encoding mechanistic relationships between variables, ODE models offer interpretability, support reasoning about underlying processes, and enable extrapolation beyond observed data\. These properties make ODEs particularly valuable in clinical settings, where understanding disease mechanisms is critical\.

However, constructing ODE models typically requires substantial expert knowledge to specify both the functional structure and parameters of the system\. In rare diseases, such knowledge is often limited, while available data are typically scarce, noisy, heterogeneous, and subject to privacy constraints\(Hilgerset al\.,[2016](https://arxiv.org/html/2607.00733#bib.bib1)\)\. Capturing this heterogeneity requires modeling parameters as distributions rather than fixed values\(Jaqaman and Danuser,[2006](https://arxiv.org/html/2607.00733#bib.bib15)\)\. These challenges make reliable model identification difficult in practice, even for approaches that rely on fitting individual\-level data\. In practice, population\-level summary statistics are often used as a pragmatic alternative, as they enable privacy\-preserving data sharing and can still provide an informative representation under noise, irregular sampling, and small cohort sizes\.

While automatic structure discovery methods such as symbolic regression\(Makke and Chawla,[2024](https://arxiv.org/html/2607.00733#bib.bib17)\)and sparse identification\(Bruntonet al\.,[2016](https://arxiv.org/html/2607.00733#bib.bib16)\)can recover governing equations from data, they rely on dense individual\-level data for model fitting and evaluation\. Recent LLM\-based ODE discovery approaches\(Holtet al\.,[2024](https://arxiv.org/html/2607.00733#bib.bib19); Duet al\.,[2024](https://arxiv.org/html/2607.00733#bib.bib20); Yanget al\.,[2026](https://arxiv.org/html/2607.00733#bib.bib21); Bideh and Gryak,[2026](https://arxiv.org/html/2607.00733#bib.bib22)\)extend this paradigm by leveraging scientific knowledge embedded in LLMs to propose candidate structures, but similarly evaluate fitness against individual\-level data\. Meanwhile, population\-level approaches such as nonlinear mixed\-effects models\(Lavielle,[2014](https://arxiv.org/html/2607.00733#bib.bib18)\)address data sparsity by pooling information across individuals, but require a pre\-specified model structure\. Together, these limitations expose a critical gap: no existing method jointly discovers ODE structure and refines parameter distributions from population\-level summary statistics alone\.

Large language models have demonstrated strong capabilities in scientific knowledge retrieval, mathematical reasoning, and code generation\(Romera\-Paredeset al\.,[2024](https://arxiv.org/html/2607.00733#bib.bib25); Shojaeeet al\.,[2024](https://arxiv.org/html/2607.00733#bib.bib5)\), as well as visual interpretation of complex data\(Achiamet al\.,[2023](https://arxiv.org/html/2607.00733#bib.bib24)\)and proposal of informative prior distributions in scientific contexts\(Riegleret al\.,[2025](https://arxiv.org/html/2607.00733#bib.bib4)\)\. These capabilities map naturally onto the core challenges of ODE discovery in data\-scarce settings\. Scientific knowledge retrieval and code generation enable the proposal of mechanistically grounded ODE structures without explicit expert input\. Visual reasoning allows the model to interpret patterns from summary statistics presented as visual summaries\. The ability to propose informative prior distributions motivates modeling parameters as distributions to capture population\-level heterogeneity\. Finally, the agentic paradigm enables LLMs to interleave reasoning and tool use, iteratively diagnosing and refining outputs based on feedback\(Yaoet al\.,[2022](https://arxiv.org/html/2607.00733#bib.bib23)\)\.

Building on these insights, we propose AgentODE \(Figure[1](https://arxiv.org/html/2607.00733#S1.F1)\), an end\-to\-end framework that jointly discovers ODE structures and refines parameter distributions from population\-level summary statistics\. At its core, AgentODE consists of two coupled feedback loops mediated by the LLM: an outer structure\-discovery loop that proposes and scores candidate ODE structures based on discrepancies between synthetic and empirical summary statistics, and an inner parameter\-inference loop that iteratively refines parameter distributions by diagnosing and reducing these discrepancies\. The process begins with the LLM proposing candidate ODE structures\. For each structure, the parameter inference agent refines parameter distributions through an iterative diagnosis–update loop\. At each iteration, the current distributions are used to simulate synthetic data, from which summary statistics are computed and compared with empirical summary statistics to identify discrepancies\. The agent then reasons over these discrepancies together with memory of previous high\-performing iterations to update the parameter distributions\. An experience buffer with dynamic memory stores discovered results and provides examples to guide subsequent ODE structure proposals\.

We evaluate AgentODE on three synthetic benchmark problems spanning materials science, pharmacokinetics, and molecular biology, as well as two clinical datasets: acute kidney injury \(AKI\) and recessive dystrophic epidermolysis bullosa \(RDEB\), a rare skin disorder with only 231 observations across 46 patients\. Using GPT\-5\.2 as the backbone LLM, AgentODE achieves performance comparable to baselines with full trajectory access, while operating on population\-level summary statistics\. On RDEB, we show that reasoning from population\-level summary statistics provides a more robust signal for structure discovery: by condensing population\-level patterns while mitigating individual\-level noise, summary statistics guide the LLM toward the true underlying mechanism, whereas baselines with access to individual\-level data, despite achieving lower prediction error, discover mechanistically implausible or degenerate structures\. This demonstrates that in sparse and noisy clinical settings, predictive performance alone is an insufficient criterion for evaluating mechanistic model quality\.Our contributionsare as follows:

- •We propose AgentODE, the first framework that jointly discovers ODE structures and refines parameter distributions from population\-level summary statistics\.
- •We demonstrate that, in sparse and noisy clinical settings, summary statistics can serve as a robust and informative data representation to guide mechanistic model inference, enabling recovery of underlying dynamical behavior while mitigating overfitting to individual\-level noise\.
- •We show that the jointly inferred structure and parameter distributions can be used to generate synthetic patient populations that capture observed heterogeneity\.

![Refer to caption](https://arxiv.org/html/2607.00733v1/x1.png)Figure 1:Overview ofAgentODE\. Inputs consist of a problem specification and empirical summaries\. An LLM proposes candidate ODE structures\. For each structure, the agent performs aninitial inferenceto obtain starting parameter distributions, from which synthetic data are generated, evaluated using logSL, and used to produce comparative visual summaries against empirical data\. The agent then refines the parameters through adiagnosis–updateloop\. Indiagnosis, discrepancies are identified from comparison visuals and quantified using statistical tools to produce a diagnostic report\. Inupdate, the agent reads this report and learns from past parameter distributions and corresponding logSL information to update the parameter distributions\. Anexperience bufferstores structures and parameter distributions to guide future structure proposals\. AgentODE outputs the best ODE modelf∗f^\{\*\}and inferred parameter distributionsp∗​\(θ\)p^\{\*\}\(\\theta\)\.
## 2Related Work

Agentic reasoning with LLMs\.The ReAct framework\(Yaoet al\.,[2022](https://arxiv.org/html/2607.00733#bib.bib23)\)established the paradigm of interleaving reasoning and tool use in LLM agents\. Reflexion\(Shinnet al\.,[2023](https://arxiv.org/html/2607.00733#bib.bib26)\)extended this idea by incorporating reflection over past failures stored in an episodic memory buffer, allowing agents to improve decisions across trials\. Related approaches have further explored structured reasoning and planning in LLMs\(Yaoet al\.,[2023](https://arxiv.org/html/2607.00733#bib.bib29); Wanget al\.,[2023](https://arxiv.org/html/2607.00733#bib.bib30)\)\. LLMs have also demonstrated strong capabilities in scientific tasks such as knowledge retrieval\(Romera\-Paredeset al\.,[2024](https://arxiv.org/html/2607.00733#bib.bib25); Shojaeeet al\.,[2024](https://arxiv.org/html/2607.00733#bib.bib5)\), visual interpretation of complex data\(Achiamet al\.,[2023](https://arxiv.org/html/2607.00733#bib.bib24)\), and the proposal of informative prior distributions\(Riegleret al\.,[2025](https://arxiv.org/html/2607.00733#bib.bib4)\)\. Building on these advances, AgentODE applies agentic reasoning to mechanistic model discovery, where iterative feedback from summary statistics guides both structure discovery and parameter inference\.

Equation and ODE structure discovery\.Classical methods including symbolic regression\(Makke and Chawla,[2024](https://arxiv.org/html/2607.00733#bib.bib17)\)and sparse identification\(Bruntonet al\.,[2016](https://arxiv.org/html/2607.00733#bib.bib16)\)recover governing equations from data but require individual\-level data and rely on predefined operator sets or function libraries\. Related advances in scientific machine learning, such as neural ODEs\(Chenet al\.,[2018](https://arxiv.org/html/2607.00733#bib.bib32)\)and universal differential equations\(Rackauckaset al\.,[2020](https://arxiv.org/html/2607.00733#bib.bib33)\), enable flexible modeling of dynamical systems but rely on uninterpretable neural components rather than explicit mechanistic structures\. LLM\-based approaches such as LLM\-SR\(Shojaeeet al\.,[2024](https://arxiv.org/html/2607.00733#bib.bib5)\)and LaSR\(Grayeliet al\.,[2024](https://arxiv.org/html/2607.00733#bib.bib34)\)partially address this gap by leveraging scientific knowledge embedded in LLMs to guide equation discovery through search or reusable concept libraries\. Recent approaches including LLM4ED\(Duet al\.,[2024](https://arxiv.org/html/2607.00733#bib.bib20)\), D3\(Holtet al\.,[2024](https://arxiv.org/html/2607.00733#bib.bib19)\), KeplerAgent\(Yanget al\.,[2026](https://arxiv.org/html/2607.00733#bib.bib21)\)leverage LLMs to propose ODE structures but still depend on individual\-level data for evaluation\. Population\-level approaches such as nonlinear mixed\-effects models\(Lavielle,[2014](https://arxiv.org/html/2607.00733#bib.bib18)\)pool information across individuals but require a predefined ODE structure\. To our knowledge, AgentODE is the first framework to jointly discover ODE structure and refine parameter distributions from population\-level summary statistics\.

Simulation\-based inference\.Simulation\-based inference encompasses likelihood\-free parameter estimation methods\(Cranmeret al\.,[2020](https://arxiv.org/html/2607.00733#bib.bib3)\)\. Classical approaches such as approximate Bayesian computation\(Beaumont,[2019](https://arxiv.org/html/2607.00733#bib.bib35)\)compare empirical and synthetic data through summary statistics, while more recent neural methods such as sequential neural likelihood\(Papamakarioset al\.,[2019](https://arxiv.org/html/2607.00733#bib.bib36)\)and neural posterior estimation\(Greenberget al\.,[2019](https://arxiv.org/html/2607.00733#bib.bib38)\)learn posterior distributions using neural density estimators\. These approaches can require large simulation budgets and can produce less interpretable outputs\. The synthetic likelihood\(Wood,[2010](https://arxiv.org/html/2607.00733#bib.bib2)\)defines a likelihood over summary statistics, enabling parameter inference by comparing synthetic and empirical summary statistics\. This aligns with our setting, where parameter inference relies on summary statistics\.

Table[1](https://arxiv.org/html/2607.00733#S2.T1)compares existing methods and AgentODE across key capabilities\. LLM\-SR†denotes our adaptation for ODE discovery\. AgentODE is the only method that jointly discovers ODE structure, infers parameter distributionsp​\(θ\)p\(\\theta\), and operates on population\-level summary statistics\.

Table 1:Method comparison across key capabilities\.†Adaptations for ODE discovery\.

## 3AgentODE

AgentODE is designed for clinical settings such as rare diseases, operating on population\-level summary statistics rather than directly using individual\-level data\. The framework consists of two tightly coupled components: an LLM\-guided ODE structure discovery module \(the outer loop; Section[3\.4](https://arxiv.org/html/2607.00733#S3.SS4)\) and a tool\-augmented parameter inference agent \(the inner loop; Section[3\.5](https://arxiv.org/html/2607.00733#S3.SS5)\)\. The structure discovery module proposes candidate ODE structures, which are evaluated by the parameter inference agent through iterative refinement of parameter distributions\. This evaluation is based on discrepancies between synthetic and empirical summary statistics, which guide subsequent structure proposals, forming a closed\-loop process that jointly improves both structure discovery and parameter inference\. We first formally introduce the problem setting\.

### 3\.1Problem Formulation

Our goal is to infer an ODE structuref∗f^\{\*\}and parameter distributionp∗​\(θ\)p^\{\*\}\(\\theta\)such that synthetic data generated by the model produce summary statistics consistent with the empirical data, thereby capturing the underlying dynamical patterns of the system\.

We consider ODEs of the formd​Xd​t=f​\(X,θ\)\\frac\{dX\}\{dt\}=f\(X,\\theta\), whereX∈ℝdX\\in\\mathbb\{R\}^\{d\}denotes the state variables andθ∈ℝm\\theta\\in\\mathbb\{R\}^\{m\}the model parameters\. To capture inter\-individual heterogeneity, we assume independent LogNormal distributions for parameters,θj∼LogNormal​\(μj,σj2\)\\theta\_\{j\}\\sim\\mathrm\{LogNormal\}\(\\mu\_\{j\},\\sigma\_\{j\}^\{2\}\)\. Given empirical data𝒳=\{xi​\(t\)\}i=1N\\mathcal\{X\}=\\\{x\_\{i\}\(t\)\\\}\_\{i=1\}^\{N\}consisting ofNNindividual trajectories, we represent the data using population\-level summary statisticssemp=ϕ​\(𝒳\)∈ℝqs\_\{\\mathrm\{emp\}\}=\\phi\(\\mathcal\{X\}\)\\in\\mathbb\{R\}^\{q\}, whereϕ\\phidenotes a summary\-statistic mapping andqqis the number of summary statistics\. For a candidate structureffand parameter distributionp​\(θ\)p\(\\theta\), synthetic data𝒳~​\(f,p​\(θ\)\)\\tilde\{\\mathcal\{X\}\}\(f,p\(\\theta\)\)are generated, from which summary statisticsssyn=ϕ​\(𝒳~​\(f,p​\(θ\)\)\)s\_\{\\mathrm\{syn\}\}=\\phi\(\\tilde\{\\mathcal\{X\}\}\(f,p\(\\theta\)\)\)are computed\. The objective is to inferf∗f^\{\*\}andp∗​\(θ\)p^\{\*\}\(\\theta\)by minimizing a discrepancy between empirical and synthetic summary statistics,\(f∗,p∗​\(θ\)\)=arg⁡minf,p​\(θ\)⁡𝒟​\(ssyn,semp\)\(f^\{\*\},p^\{\*\}\(\\theta\)\)=\\arg\\min\_\{f,\\,p\(\\theta\)\}\\mathcal\{D\}\\big\(s\_\{\\mathrm\{syn\}\},s\_\{\\mathrm\{emp\}\}\\big\), where𝒟​\(⋅,⋅\)\\mathcal\{D\}\(\\cdot,\\cdot\)denotes a discrepancy function defined over summary statistics, as detailed in Sections[3\.2](https://arxiv.org/html/2607.00733#S3.SS2)and[3\.3](https://arxiv.org/html/2607.00733#S3.SS3)\.

### 3\.2Summary Statistics

Summary statistics support both structure discovery and parameter inference\. For structure discovery, we use a broad set of statistics capturing distributions, autocorrelation, stationarity, entropy, scaling, trend, volatility, and cross\-variable dependencies\(Fulcher,[2018](https://arxiv.org/html/2607.00733#bib.bib6); Fulcher and Jones,[2014](https://arxiv.org/html/2607.00733#bib.bib7); Bandt and Pompe,[2002](https://arxiv.org/html/2607.00733#bib.bib8); Kantelhardtet al\.,[2001](https://arxiv.org/html/2607.00733#bib.bib9); Hyndman and Athanasopoulos,[2018](https://arxiv.org/html/2607.00733#bib.bib10)\), enabling informative comparison between candidate models\. For parameter inference, we use a focused subset of the summary statistics used for structure comparison to provide stable signals for likelihood estimation\. To better capture the dynamics of sparse and noisy clinical data, we prioritize distributional shape, temporal autocorrelation, population\-level trends, and cross\-variable dependencies, which together reflect the most informative aspects of population dynamics while ensuring consistency between parameter inference and structure comparison \(complete set of statistics is provided in Appendix[A](https://arxiv.org/html/2607.00733#A1)\)\.

### 3\.3Scoring

AgentODE employs two complementary metrics: log synthetic likelihood \(logSL\) for parameter inference within an ODE structure, and mean normalized summary discrepancy \(MNSD\) for comparison across ODE structures\. LogSL evaluates parameter fit within a structure but is not comparable across structures, as each structure has its own parameterization and induces a different distribution over summary statistics\. We therefore introduce MNSD for consistent comparison across structures\.

Scoring within a structure \(logSL\)\.Given a candidate structureffand parameter distributionp​\(θ\)p\(\\theta\)parameterized byψ=\{\(μj,σj\)\}j=1m\\psi=\\\{\(\\mu\_\{j\},\\sigma\_\{j\}\)\\\}\_\{j=1\}^\{m\}, we performNsimN\_\{\\mathrm\{sim\}\}simulations by drawing independent parameter vectors\{θ\(n\)\}n=1Nsim\\\{\\theta^\{\(n\)\}\\\}\_\{n=1\}^\{N\_\{\\mathrm\{sim\}\}\}fromp​\(θ\)p\(\\theta\)\. Each simulation generates a population of synthetic trajectories𝒳~\(n\)​\(f,θ\(n\)\)\\tilde\{\\mathcal\{X\}\}^\{\(n\)\}\(f,\\theta^\{\(n\)\}\), from which summary statisticsssynSL,\(n\)=ϕSL​\(𝒳~\(n\)​\(f,θ\(n\)\)\)s\_\{\\mathrm\{syn\}\}^\{\\mathrm\{SL\},\(n\)\}=\\phi\_\{\\mathrm\{SL\}\}\(\\tilde\{\\mathcal\{X\}\}^\{\(n\)\}\(f,\\theta^\{\(n\)\}\)\)are computed\. Following\(Wood,[2010](https://arxiv.org/html/2607.00733#bib.bib2)\), we evaluate model fit using the log synthetic likelihoodlog⁡SL​\(ψ\)=−12​\[log⁡\|Σ^ψ\|\+\(sempSL−μ^ψ\)⊤​Σ^ψ−1​\(sempSL−μ^ψ\)\]\\log\\mathrm\{SL\}\(\\psi\)=\-\\tfrac\{1\}\{2\}\\left\[\\log\|\\hat\{\\Sigma\}\_\{\\psi\}\|\+\(s\_\{\\mathrm\{emp\}\}^\{\\mathrm\{SL\}\}\-\\hat\{\\mu\}\_\{\\psi\}\)^\{\\top\}\\hat\{\\Sigma\}\_\{\\psi\}^\{\-1\}\(s\_\{\\mathrm\{emp\}\}^\{\\mathrm\{SL\}\}\-\\hat\{\\mu\}\_\{\\psi\}\)\\right\], whereμ^ψ\\hat\{\\mu\}\_\{\\psi\}andΣ^ψ\\hat\{\\Sigma\}\_\{\\psi\}denote the mean and covariance of\{ssynSL,\(n\)\}n=1Nsim\\\{s\_\{\\mathrm\{syn\}\}^\{\\mathrm\{SL\},\(n\)\}\\\}\_\{n=1\}^\{N\_\{\\mathrm\{sim\}\}\}\. The logSL score is maximized during parameter inference\.

Scoring across structures \(MNSD\)\.To compare candidate ODE structures, we quantify discrepancies between empirical and synthetic summary statistics using a normalized absolute difference, which we term the mean normalized summary discrepancy \(MNSD\):MNSD=1\|ℐ\|​∑i∈ℐ\|semp,iMNSD−ssyn,iMNSD\|IQRiemp\.\\mathrm\{MNSD\}=\\frac\{1\}\{\|\\mathcal\{I\}\|\}\\sum\_\{i\\in\\mathcal\{I\}\}\\frac\{\\left\|s^\{\\mathrm\{MNSD\}\}\_\{\\mathrm\{emp\},i\}\-s^\{\\mathrm\{MNSD\}\}\_\{\\mathrm\{syn\},i\}\\right\|\}\{\\mathrm\{IQR\}^\{\\mathrm\{emp\}\}\_\{i\}\}\.

Here,semp,iMNSDs^\{\\mathrm\{MNSD\}\}\_\{\\mathrm\{emp\},i\}denotes theii\-th empirical summary statistic, andssyn,iMNSDs^\{\\mathrm\{MNSD\}\}\_\{\\mathrm\{syn\},i\}is computed from synthetic data generated using the inferred parameter distribution for the candidate structure\.IQRiemp\\mathrm\{IQR\}^\{\\mathrm\{emp\}\}\_\{i\}is the empirical interquartile range of that statistic\(Law,[1986](https://arxiv.org/html/2607.00733#bib.bib28)\)\. The index setℐ\\mathcal\{I\}indexes the summary statistics used for structure comparison\. Lower MNSD indicates better agreement between empirical and synthetic summary statistics\.

### 3\.4ODE Structure Discovery

Base framework: LLM\-SR\.We build on LLM\-SR\(Shojaeeet al\.,[2024](https://arxiv.org/html/2607.00733#bib.bib5)\), a symbolic regression framework that iteratively generates, evaluates, and refines candidate symbolic equations using LLMs\. At each iteration, the LLM proposes equation hypotheses conditioned on the problem description, which are then evaluated based on their fit to individual\-level data\. The candidate equations and their corresponding fit scores are stored in an experience buffer, which serves as a dynamic memory to guide subsequent proposals\. The buffer adopts an island\-based model\(Tanese,[1989](https://arxiv.org/html/2607.00733#bib.bib27)\)to maintain diversity and avoid local optima\.

Extension to ODE discovery\.We extend this paradigm to dynamical systems by replacing symbolic regression with ODE structure discovery\. Unlike LLM\-SR, which directly evaluates candidate equations using individual\-level data, our framework evaluates ODE structures via a parameter inference agent[3\.5](https://arxiv.org/html/2607.00733#S3.SS5)\. For each candidate structure, the agent refines parameter distributions using logSL, after which the structure is scored using MNSD\.

### 3\.5Parameter Inference Agent

Given a candidate ODE structureff, the agent operates in two phases: an initial inference producing a starting parameter distribution, followed by an iterative diagnosis–update loop that refines the parameter distributionp​\(θ\)p\(\\theta\), guided by logSL measuring agreement betweenssynSLs\_\{\\mathrm\{syn\}\}^\{\\mathrm\{SL\}\}andsempSLs\_\{\\mathrm\{emp\}\}^\{\\mathrm\{SL\}\}\.

Visual summaries\.We construct visual summaries from population\-level data in a structured and interpretable form to represent temporal dynamics, heterogeneity, and inter\-variable relationships\. These include \(i\) mean trajectories with 95% confidence intervals over time, \(ii\) mean trajectories with 95% confidence intervals stratified by initial\-condition quartiles, and \(iii\) heatmaps of Spearman correlations of first differences between variables\. These summaries are used as empirical inputs for initial parameter inference\. During diagnosis, empirical and synthetic versions of these summaries are compared to identify discrepancies that guide parameter updates\.

LLM Initial InferenceGiven the problem specification, candidate ODE structureff, and empirical visual summaries, the LLM proposes an initial parameter distributionp\(0\)​\(θ\)p^\{\(0\)\}\(\\theta\)\. Details in Appendix[G\.2](https://arxiv.org/html/2607.00733#A7.SS2)\.

Diagnosis\.The agent diagnoses the current best\-performing parameter distribution, as measured by logSL\. It first reasons from the problem context and the ODE structureffto infer expected system dynamics, including trajectory behavior, heterogeneity, and relationships between variables\. These expectations are grounded by comparing empirical and synthetic visual summaries to identify discrepancies in the simulated dynamics\. Beyond the most apparent mismatches, the agent reasons about their underlying causes in the context offfand the problem domain, linking them to parameter roles and the dynamical assumptions implied by the current distribution\. The same summary statistics set used for structure scoring are provided as tools to the agent, which selectively evaluates them to quantify discrepancies and support its interpretation\. The results are compiled into a structured diagnostic report that characterizes each failure mode by its severity, affected variables, and supporting statistical evidence\. Details in Appendix[G\.3](https://arxiv.org/html/2607.00733#A7.SS3)\.

Update\.Given the diagnostic report and past parameter distributions with high logSL scores, the agent updates the parameter distribution\. First, it examines the diagnosed failure modes and reasons in the context of the ODE structureffand the problem domain to identify which parameters govern the affected dynamics and what directional changes may be required\. It then reviews past parameter settings and their associated logSL outcomes to understand how previous configurations influenced model fit, informing the next update and avoiding repeating ineffective configurations\.

## 4Experiments and Results

### 4\.1Benchmarks and Clinical Datasets

We evaluate AgentODE on three synthetic benchmarks and two clinical datasets, focusing on ODE structure recovery and inferred parameter distributions\.

Benchmarks\.We evaluate AgentODE on three benchmark problems spanning materials science, pharmacokinetics, and molecular biology\. To ensure evaluation beyond memorization, we assess recitation risk through contamination analysis \(Table[3](https://arxiv.org/html/2607.00733#A4.T3)in Appendix[D](https://arxiv.org/html/2607.00733#A4)\)\. TheApoptosissystem models programmed cell death dynamics across three coupled molecular species in molecular biology\(Bideh and Gryak,[2026](https://arxiv.org/html/2607.00733#bib.bib22)\)\. ThePolymer DA Cross\-linkingsystem models Diels–Alder reaction kinetics across four coupled chemical species in materials science\(Bideh and Gryak,[2026](https://arxiv.org/html/2607.00733#bib.bib22)\)\. ThePKPD\-Immunesystem couples two\-compartment pharmacokinetics with immune and effector cell dynamics across four state variables\. AgentODE is designed to operate over heterogeneous populations, requiring that benchmark systems accommodate individual\-level variability\. We therefore sample initial conditions from domain\-grounded distributions for all systems\. For parameters, we introduce inter\-individual variability where meaningful: thePKPD\-Immunesystem naturally accommodates inter\-individual variability in pharmacokinetics and immune response, and we assign log\-normal parameter distributions accordingly\. For the remaining systems, parameters are fixed following prior work\(Bideh and Gryak,[2026](https://arxiv.org/html/2607.00733#bib.bib22)\)\. Details in Appendix[D](https://arxiv.org/html/2607.00733#A4)\.

Clinical Datasets\.We evaluate AgentODE on two clinical datasets\. The first is an acute kidney injury \(AKI\) cohort extracted from MIMIC\-IV 3\.1\(Johnsonet al\.,[2023](https://arxiv.org/html/2607.00733#bib.bib39)\), comprising 353 admissions \(trajectories\) from patients aged 50–64 years, with irregularly spaced measurements of creatinine, blood urea nitrogen \(BUN\), and potassium over the first 168 hours after admission\. Measurements are irregularly spaced with variable numbers of observations per trajectory\. The second dataset comprises 46 patients with severe recessive dystrophic epidermolysis bullosa \(RDEB\)\(Bardhanet al\.,[2020](https://arxiv.org/html/2607.00733#bib.bib42); Reimeret al\.,[2020](https://arxiv.org/html/2607.00733#bib.bib40); Karakioulakiet al\.,[2026](https://arxiv.org/html/2607.00733#bib.bib41)\), a rare genetic skin disease, with 231 total observations across up to 144 months since first hospital visit, a mean of 5 observations per patient, and 16 patients with only 2 time points\. We use measurements of serum albumin, log\-transformed C\-reactive protein \(log10\(CRP\+\+1\)\), hemoglobin, and BMI relative to age\- and sex\-matched WHO reference \(BMIrel\)\.

### 4\.2Experimental Setup

We compare AgentODE against state\-of\-the\-art structure discovery baselines spanning sparse regression \(SINDy\(Bruntonet al\.,[2016](https://arxiv.org/html/2607.00733#bib.bib16)\)\), evolutionary symbolic regression \(PySR\(Cranmer,[2023](https://arxiv.org/html/2607.00733#bib.bib44)\)\), and LLM\-guided symbolic regression \(LLM\-SR\(Shojaeeet al\.,[2024](https://arxiv.org/html/2607.00733#bib.bib5)\)\), as well asNeural ODE\(Chenet al\.,[2018](https://arxiv.org/html/2607.00733#bib.bib32)\)as a black\-box baseline\. For PySR, we fit one independent symbolic regression model per state variable and assemble them into a system of equations\. LLM\-SR†denotes our adaptation of LLM\-SR to ODE discovery, where the LLM proposes ODE structure candidates governing all state variables jointly, and parameters are subsequently fitted independently for each trajectory using BFGS, providing a reference with full access to individual\-level data\. SINDy and PySR are not evaluated on clinical datasets, as they require preprocessing steps such as interpolation and derivative estimation that introduce assumptions incompatible with our setting\. We additionally include two summary\-level ablations:LLM\-based parameter initialization \(k=5k=5\), which queries the LLM 5 times given the ODE structure and problem context to propose parameter distributions and selects the one with the highest logSL for evaluation, andAgentODE w/o iterative refinement, which removes the diagnosis–update loop from AgentODE\. For parameter inference, AgentODE uses a bounded iteration budget with early stopping \(10 iterations per structure with patience 5\)\. LLM\-SR†is evaluated with over 1000 sampled candidates\. AgentODE and its ablations use approximately 500 samples per experiment\. All other baselines follow their original setups or are tuned to achieve their best performance\. Implementation details in Appendix[C](https://arxiv.org/html/2607.00733#A3)\.

### 4\.3Quantitative Results

PerformanceTable[2](https://arxiv.org/html/2607.00733#S4.T2)reports RMSE across all datasets\. Among summary\-level methods, AgentODE achieves the best performance on four of five datasets\. On benchmarks, AgentODE performs competitively with full\-trajectory methods, particularly on PKPD, where it achieves an RMSE of 0\.064 compared to 0\.685 for Neural ODE\. On clinical datasets, AgentODE substantially outperforms Neural ODE on AKI \(2\.56 vs 6\.162\) and achieves comparable performance to LLM\-SR†\. However, in sparse and noisy clinical settings, predictive error alone is insufficient for evaluating mechanistic models: lower RMSE may reflect overfitting to noisy trajectories rather than recovery of clinically meaningful dynamics\. Across datasets, AgentODE consistently recovers more compact and mechanistically meaningful ODE structures, despite operating solely on summary\-level statistics\. Qualitative comparisons of the discovered systems are provided in Appendix[F](https://arxiv.org/html/2607.00733#A6)\.

Table 2:RMSE: averaged across state variables \(lower is better\)\.Boldindicates best performance within each group\.N/A: not applicable due to data requirements\.Invalid: proposed parameter distributions yielded unstable \(e\.g\., NaN/Inf\) or physically unrealistic synthetic trajectories\. For RDEB, results are reported in\-sample due to limited cohort size\. LLM backbone:GPT\-5\.2\.BenchmarksClinicalMethodApoptosisPolymerPKPDAKIRDEBFull trajectory accessSINDy0\.0040\.5431\.251N/AN/APySR0\.0450\.1960\.981N/AN/ALLM\-SR†0\.0240\.5080\.0152\.5830\.532Neural ODE0\.1340\.2570\.6856\.1620\.528Summary\-level onlyLLM\-based parameter initialization \(k=5k=5\)InvalidInvalid0\.53InvalidInvalidAgentODE w/o iterative refinement0\.1111\.810\.2593\.8030\.628AgentODE0\.0790\.6890\.0642\.560\.644EfficiencyFigure[2](https://arxiv.org/html/2607.00733#S4.F2)shows the best score trajectory over structure discovery iterations for AgentODE, AgentODE w/o iterative refinement, and LLM\-based parameter initialization \(k=5k=5\) across all datasets\. AgentODE exhibits more consistent improvement over the search, with more frequent score improvements and fewer extended plateaus compared to both ablations, indicating that iterative refinement plays a critical role in sustaining productive exploration\. The gap between methods is most pronounced on clinical datasets, where AgentODE achieves substantially higher final scores within the same iteration budget\.

![Refer to caption](https://arxiv.org/html/2607.00733v1/x2.png)Figure 2:Ablation study comparing AgentODE, AgentODE w/o iterative refinement, and LLM\-based parameter initialization \(k=5k=5\) across five datasets\. Each panel shows the best score \(−\-MNSD, higher is better\) over the number of candidate ODE structures evaluated\. AgentODE consistently achieves higher scores and improves more efficiently than both ablations across all datasets\.Analysis\.The ablation results reveal two complementary roles of the framework design\. First, visual summaries provide an effective inductive signal for inferring valid parameter distributions, as evidenced by the Invalid results \(4/5\) for LLM\-based parameter initialization \(k=5\) across datasets\. This effect is most pronounced for Apoptosis and Polymer, where parameters are fixed and must be inferred from data patterns, and in clinical settings where visual summaries provide a direct representation of temporal patterns, without which the model produces invalid outcomes on AKI and RDEB\. For PKPD, where parameters follow distributions, embedded domain knowledge enables reasonable initialization, though visual summaries still lead to substantial improvement \(0\.064 vs 0\.259\)\. Second, diagnosis–update loop consistently improves final performance compared to AgentODE w/o iterative refinement, achieving substantially lower RMSE on Polymer \(0\.689 vs 1\.81\), PKPD \(0\.064 vs 0\.259\), and AKI \(2\.56 vs 3\.803\)\. Figure[2](https://arxiv.org/html/2607.00733#S4.F2)further shows that refinement leads to more sustained improvement during structure discovery\. Finally, benchmark problems comprise clean, fully observed synthetic trajectories, where full trajectory access provides a natural advantage, whereas clinical datasets are scarce, noisy, and heterogeneous, precisely where summary\-level inference is particularly valuable\. The RDEB case study illustrates this concretely: AgentODE w/o iterative refinement and LLM\-SR†achieve lower RMSE yet recover mechanistically implausible or degenerate ODE structures \(causal comparison in Figure[7](https://arxiv.org/html/2607.00733#A6.F7)\), suggesting that methods fitting directly to noisy individual\-level data overfit to noise rather than recovering the true underlying mechanism\. Since high\-scoring structures guide future proposals through the experience buffer, such degenerate structures can mislead the entire discovery process, as evidenced by the final structures discovered by these baselines \( Appendix[F\.5](https://arxiv.org/html/2607.00733#A6.SS5)\)\. Reasoning from population\-level summary statistics provides a more robust signal that mitigates this risk\.

### 4\.4Qualitative Results

Diagnosis–Update LoopFigure[3](https://arxiv.org/html/2607.00733#S4.F3)illustrates the diagnosis–update loop for a candidate ODE structure on the Polymer benchmark\. The logSL score improves substantially from iteration 1 to the best iteration 4\. This improvement is reflected in the visual summaries: at iteration 1, the synthetic correlation structure poorly matches the empirical heatmap, and the mean trajectories with 95% CI exhibit clear misalignment\. By iteration 4, the synthetic correlation structure closely matches the empirical one, and the mean trajectories with 95% CI align better with the empirical data\.

Analysis\.The diagnosis–update loop improves parameter distributions, as reflected by consistent logSL improvements in early iterations \(Appendix[G\.3](https://arxiv.org/html/2607.00733#A7.SS3)\)\. However, the non\-monotonic behavior, where the score slightly decreases after the best iteration, highlights an inherent limitation: parameter distributions are updated through LLM\-based reasoning over few iterations, rather than gradient\-based optimization, and thus monotonic convergence is not guaranteed\. This is further compounded by parameter non\-identifiability\(Raueet al\.,[2009](https://arxiv.org/html/2607.00733#bib.bib49)\), where multiple parameter configurations may yield similarly plausible summary statistics\. Nevertheless, for a substantial subset of parameters, the inferred distributions fall within plausible ranges and are more concentrated than empirically fitted distributions obtained via trajectory\-level optimization \(Figure[8](https://arxiv.org/html/2607.00733#A6.F8)\)\. This behavior is consistent with LLM\-based inference concentrating probability mass around plausible values when reasoning from population\-level summary statistics, as further supported by correlation heatmaps \(Figure[9](https://arxiv.org/html/2607.00733#A6.F9)\)\.

![Refer to caption](https://arxiv.org/html/2607.00733v1/x3.png)Figure 3:Diagnosis\-update iterationsfor an ODE structure on polymer\.Left:logSL score across iterations, with the best iteration marked\.Right:Spearman correlation heatmaps and mean trajectories±\\pm95% CI comparing empirical and synthetic data at iteration 1 and the best iteration\.

## 5Conclusion

We introduced AgentODE, an end\-to\-end framework that jointly discovers ODE structures and infers parameter distributions from population\-level summary statistics\. Across benchmarks and clinical datasets, AgentODE achieves performance comparable to baselines that rely on individual\-level data\. The rare disease RDEB case study demonstrates that reasoning from summary statistics provides a robust signal for structure discovery, mitigating overfitting to noisy clinical observations\. AgentODE recovers a mechanistically consistent structure aligned with an independently developed expert model, whereas baselines with full data access produce structurally less coherent models that deviate from clinically meaningful relationships\. Crucially, operating on population\-level summary statistics provides a principled privacy\-preserving way to mechanistic modeling for rare disease research\.

Limitations and Broader Impact\.AgentODE’s performance relies on knowledge embedded in the LLM and decreases with increasing ODE complexity\. It is not intended to replace numerical optimization with full trajectories; parameter identifiability and distribution refinement under data scarcity remain fundamental challenges\. Rather, AgentODE provides a practical approach for settings where data are limited, noisy, heterogeneous, or privacy\-constrained\. AgentODE has direct applications in biomedical research for under\-studied diseases where privacy constraints limit access to individual\-level data\. Recovered structures and distributions should be interpreted as functionally consistent with observed population behavior rather than ground truth, and validated by domain experts before clinical use\.

Future Work and Outlook\.AgentODE opens several promising directions\. On the methodological side, Meta\-Harness–style approaches\(Leeet al\.,[2026](https://arxiv.org/html/2607.00733#bib.bib47)\)could enable automatic optimization of the agentic pipeline, improving parameter distribution refinement without manual design\. Extending the framework to a society of thought\(Evanset al\.,[2026](https://arxiv.org/html/2607.00733#bib.bib48)\), where multiple specialized agents jointly reason over structure and parameters, could further strengthen the diagnosis–update loop\. We aim to bring AgentODE into practice to help clinicians understand and cure rare diseases in settings with data scarcity and privacy constraints\.

## Acknowledgments and Disclosure of Funding

Funded by the Deutsche Forschungsgemeinschaft \(DFG, German Research Foundation\) – Project\-ID 499552394 – SFB 1597\.

## References

- J\. Achiam, S\. Adler, S\. Agarwal, L\. Ahmad, I\. Akkaya, F\. L\. Aleman, D\. Almeida, J\. Altenschmidt, S\. Altman, S\. Anadkat,et al\.\(2023\)Gpt\-4 technical report\.arXiv preprint arXiv:2303\.08774\.Cited by:[§1](https://arxiv.org/html/2607.00733#S1.p4.1),[§2](https://arxiv.org/html/2607.00733#S2.p1.1)\.
- Permutation entropy: a natural complexity measure for time series\.Physical review letters88\(17\),pp\. 174102\.Cited by:[§3\.2](https://arxiv.org/html/2607.00733#S3.SS2.p1.1)\.
- A\. Bardhan, L\. Bruckner\-Tuderman, I\. L\. Chapple, J\. Fine, N\. Harper, C\. Has, T\. M\. Magin, M\. P\. Marinkovich, J\. F\. Marshall, J\. A\. McGrath,et al\.\(2020\)Epidermolysis bullosa\.Nature Reviews Disease Primers6\(1\),pp\. 78\.Cited by:[§E\.2](https://arxiv.org/html/2607.00733#A5.SS2.p1.1),[§4\.1](https://arxiv.org/html/2607.00733#S4.SS1.p3.2)\.
- M\. A\. Beaumont \(2019\)Approximate bayesian computation\.Annual review of statistics and its application6\(1\),pp\. 379–403\.Cited by:[§2](https://arxiv.org/html/2607.00733#S2.p3.1)\.
- A\. Z\. Bideh and J\. Gryak \(2026\)LLM\-ode: data\-driven discovery of dynamical systems with large language models\.arXiv preprint arXiv:2603\.20910\.Cited by:[§D\.1](https://arxiv.org/html/2607.00733#A4.SS1.p1.10),[§D\.2](https://arxiv.org/html/2607.00733#A4.SS2.p1.4),[§D\.2](https://arxiv.org/html/2607.00733#A4.SS2.p4.6),[§1](https://arxiv.org/html/2607.00733#S1.p3.1),[§4\.1](https://arxiv.org/html/2607.00733#S4.SS1.p2.1)\.
- S\. L\. Brunton, J\. L\. Proctor, and J\. N\. Kutz \(2016\)Discovering governing equations from data by sparse identification of nonlinear dynamical systems\.Proceedings of the national academy of sciences113\(15\),pp\. 3932–3937\.Cited by:[Appendix B](https://arxiv.org/html/2607.00733#A2.p2.2),[§1](https://arxiv.org/html/2607.00733#S1.p3.1),[§2](https://arxiv.org/html/2607.00733#S2.p2.1),[§4\.2](https://arxiv.org/html/2607.00733#S4.SS2.p1.3)\.
- R\. T\. Chen, Y\. Rubanova, J\. Bettencourt, and D\. K\. Duvenaud \(2018\)Neural ordinary differential equations\.Advances in neural information processing systems31\.Cited by:[Appendix B](https://arxiv.org/html/2607.00733#A2.p5.2),[§2](https://arxiv.org/html/2607.00733#S2.p2.1),[§4\.2](https://arxiv.org/html/2607.00733#S4.SS2.p1.3)\.
- K\. Cranmer, J\. Brehmer, and G\. Louppe \(2020\)The frontier of simulation\-based inference\.Proceedings of the National Academy of Sciences117\(48\),pp\. 30055–30062\.Cited by:[§2](https://arxiv.org/html/2607.00733#S2.p3.1)\.
- M\. Cranmer \(2023\)Interpretable machine learning for science with pysr and symbolicregression\. jl\.arXiv preprint arXiv:2305\.01582\.Cited by:[Appendix B](https://arxiv.org/html/2607.00733#A2.p3.2),[§4\.2](https://arxiv.org/html/2607.00733#S4.SS2.p1.3)\.
- B\. M\. De Silva, K\. Champion, M\. Quade, J\. Loiseau, J\. N\. Kutz, and S\. L\. Brunton \(2020\)Pysindy: a python package for the sparse identification of nonlinear dynamics from data\.arXiv preprint arXiv:2004\.08424\.Cited by:[Appendix B](https://arxiv.org/html/2607.00733#A2.p2.2)\.
- M\. Du, Y\. Chen, Z\. Wang, L\. Nie, and D\. Zhang \(2024\)Llm4ed: large language models for automatic equation discovery\.arXiv preprint arXiv:2405\.07761\.Cited by:[§1](https://arxiv.org/html/2607.00733#S1.p3.1),[§2](https://arxiv.org/html/2607.00733#S2.p2.1)\.
- J\. Evans, B\. Bratton, and B\. Agüera y Arcas \(2026\)Agentic ai and the next intelligence explosion\.Vol\.391,American Association for the Advancement of Science\.Cited by:[§5](https://arxiv.org/html/2607.00733#S5.p3.1)\.
- B\. D\. Fulcher and N\. S\. Jones \(2014\)Highly comparative feature\-based time\-series classification\.IEEE Transactions on Knowledge and Data Engineering26\(12\),pp\. 3026–3037\.Cited by:[§3\.2](https://arxiv.org/html/2607.00733#S3.SS2.p1.1)\.
- B\. D\. Fulcher \(2018\)Feature\-based time\-series analysis\.InFeature engineering for machine learning and data analytics,pp\. 87–116\.Cited by:[§3\.2](https://arxiv.org/html/2607.00733#S3.SS2.p1.1)\.
- A\. Grayeli, A\. Sehgal, O\. Costilla\-Reyes, M\. Cranmer, and S\. Chaudhuri \(2024\)Symbolic regression with a learned concept library\.Advances in Neural Information Processing Systems37,pp\. 44678–44709\.Cited by:[§2](https://arxiv.org/html/2607.00733#S2.p2.1)\.
- D\. Greenberg, M\. Nonnenmacher, and J\. Macke \(2019\)Automatic posterior transformation for likelihood\-free inference\.InInternational conference on machine learning,pp\. 2404–2414\.Cited by:[§2](https://arxiv.org/html/2607.00733#S2.p3.1)\.
- R\. Hilgers, F\. König, G\. Molenberghs, and S\. Senn \(2016\)Design and analysis of clinical trials for small rare disease populations\.Journal of Rare Diseases Research & Treatment1\(1\),pp\. 53–60\.Cited by:[§1](https://arxiv.org/html/2607.00733#S1.p2.1)\.
- S\. Holt, Z\. Qian, T\. Liu, J\. Weatherall, and M\. van der Schaar \(2024\)Data\-driven discovery of dynamical systems in pharmacology using large language models\.Advances in Neural Information Processing Systems37,pp\. 96325–96366\.Cited by:[§1](https://arxiv.org/html/2607.00733#S1.p3.1),[§2](https://arxiv.org/html/2607.00733#S2.p2.1)\.
- R\. J\. Hyndman and G\. Athanasopoulos \(2018\)Forecasting: principles and practice\.OTexts\.Cited by:[§3\.2](https://arxiv.org/html/2607.00733#S3.SS2.p1.1)\.
- K\. Jaqaman and G\. Danuser \(2006\)Linking data to models: data regression\.Nature Reviews Molecular Cell Biology7\(11\),pp\. 813–819\.Cited by:[§1](https://arxiv.org/html/2607.00733#S1.p2.1)\.
- A\. E\. Johnson, L\. Bulgarelli, L\. Shen, A\. Gayles, A\. Shammout, S\. Horng, T\. J\. Pollard, S\. Hao, B\. Moody, B\. Gow,et al\.\(2023\)MIMIC\-iv, a freely accessible electronic health record dataset\.Scientific data10\(1\),pp\. 1\.Cited by:[§E\.1](https://arxiv.org/html/2607.00733#A5.SS1.p1.4),[§4\.1](https://arxiv.org/html/2607.00733#S4.SS1.p3.2)\.
- J\. W\. Kantelhardt, E\. Koscielny\-Bunde, H\. H\. Rego, S\. Havlin, and A\. Bunde \(2001\)Detecting long\-range correlations with detrended fluctuation analysis\.Physica A: Statistical Mechanics and its Applications295\(3\-4\),pp\. 441–454\.Cited by:[§3\.2](https://arxiv.org/html/2607.00733#S3.SS2.p1.1)\.
- M\. Karakioulaki, H\. Yang, M\. Hess, H\. Binder, K\. Eyerich, and C\. Has \(2026\)Systemic inflammation in recessive dystrophic epidermolysis bullosa: a five\-year longitudinal study\.British Journal of Dermatology,pp\. ljag082\.Cited by:[§E\.2](https://arxiv.org/html/2607.00733#A5.SS2.p1.1),[§4\.1](https://arxiv.org/html/2607.00733#S4.SS1.p3.2)\.
- D\. P\. Kingma and J\. Ba \(2014\)Adam: a method for stochastic optimization\.arXiv preprint arXiv:1412\.6980\.Cited by:[Appendix B](https://arxiv.org/html/2607.00733#A2.p5.2)\.
- M\. Lavielle \(2014\)Mixed effects models for the population approach: models, tasks, methods and tools\.CRC press\.Cited by:[§1](https://arxiv.org/html/2607.00733#S1.p3.1),[§2](https://arxiv.org/html/2607.00733#S2.p2.1)\.
- J\. Law \(1986\)Robust statistics—the approach based on influence functions\.Wiley Online Library\.Cited by:[§3\.3](https://arxiv.org/html/2607.00733#S3.SS3.p4.5)\.
- Y\. Lee, R\. Nair, Q\. Zhang, K\. Lee, O\. Khattab, and C\. Finn \(2026\)Meta\-harness: end\-to\-end optimization of model harnesses\.arXiv preprint arXiv:2603\.28052\.Cited by:[§5](https://arxiv.org/html/2607.00733#S5.p3.1)\.
- N\. Makke and S\. Chawla \(2024\)Interpretable scientific discovery with symbolic regression: a review\.Artificial Intelligence Review57\(1\),pp\. 2\.Cited by:[§1](https://arxiv.org/html/2607.00733#S1.p3.1),[§2](https://arxiv.org/html/2607.00733#S2.p2.1)\.
- G\. Papamakarios, D\. Sterratt, and I\. Murray \(2019\)Sequential neural likelihood: fast likelihood\-free inference with autoregressive flows\.InThe 22nd international conference on artificial intelligence and statistics,pp\. 837–848\.Cited by:[§2](https://arxiv.org/html/2607.00733#S2.p3.1)\.
- C\. Rackauckas, Y\. Ma, J\. Martensen, C\. Warner, K\. Zubov, R\. Supekar, D\. Skinner, A\. Ramadhan, and A\. Edelman \(2020\)Universal differential equations for scientific machine learning\.arXiv preprint arXiv:2001\.04385\.Cited by:[§2](https://arxiv.org/html/2607.00733#S2.p2.1)\.
- A\. Raue, C\. Kreutz, T\. Maiwald, J\. Bachmann, M\. Schilling, U\. Klingmüller, and J\. Timmer \(2009\)Structural and practical identifiability analysis of partially observed dynamical models by exploiting the profile likelihood\.Bioinformatics25\(15\),pp\. 1923–1929\.Cited by:[§4\.4](https://arxiv.org/html/2607.00733#S4.SS4.p2.1)\.
- A\. Raue, M\. Schilling, J\. Bachmann, A\. Matteson, M\. Schelke, D\. Kaschek, S\. Hug, C\. Kreutz, B\. D\. Harms, F\. J\. Theis,et al\.\(2013\)Lessons learned from quantitative dynamical modeling in systems biology\.PloS one8\(9\),pp\. e74335\.Cited by:[§1](https://arxiv.org/html/2607.00733#S1.p1.1)\.
- A\. Reimer, M\. Hess, A\. Schwieger\-Briel, D\. Kiritsi, F\. Schauer, H\. Schumann, L\. Bruckner\-Tuderman, and C\. Has \(2020\)Natural history of growth and anaemia in children with epidermolysis bullosa: a retrospective cohort study\.British Journal of Dermatology182\(6\),pp\. 1437–1448\.Cited by:[§E\.2](https://arxiv.org/html/2607.00733#A5.SS2.p1.1),[§4\.1](https://arxiv.org/html/2607.00733#S4.SS1.p3.2)\.
- M\. A\. Riegler, K\. H\. Hellton, V\. Thambawita, and H\. L\. Hammer \(2025\)Using large language models to suggest informative prior distributions in bayesian regression analysis\.Scientific Reports15\(1\),pp\. 33386\.Cited by:[§1](https://arxiv.org/html/2607.00733#S1.p4.1),[§2](https://arxiv.org/html/2607.00733#S2.p1.1)\.
- B\. Romera\-Paredes, M\. Barekatain, A\. Novikov, M\. Balog, M\. P\. Kumar, E\. Dupont, F\. J\. Ruiz, J\. S\. Ellenberg, P\. Wang, O\. Fawzi,et al\.\(2024\)Mathematical discoveries from program search with large language models\.Nature625\(7995\),pp\. 468–475\.Cited by:[§1](https://arxiv.org/html/2607.00733#S1.p4.1),[§2](https://arxiv.org/html/2607.00733#S2.p1.1)\.
- N\. Shinn, F\. Cassano, A\. Gopinath, K\. Narasimhan, and S\. Yao \(2023\)Reflexion: language agents with verbal reinforcement learning\.Advances in neural information processing systems36,pp\. 8634–8652\.Cited by:[§2](https://arxiv.org/html/2607.00733#S2.p1.1)\.
- P\. Shojaee, K\. Meidani, S\. Gupta, A\. B\. Farimani, and C\. K\. Reddy \(2024\)Llm\-sr: scientific equation discovery via programming with large language models\.arXiv preprint arXiv:2404\.18400\.Cited by:[Appendix B](https://arxiv.org/html/2607.00733#A2.p4.1),[Appendix C](https://arxiv.org/html/2607.00733#A3.p2.2),[§1](https://arxiv.org/html/2607.00733#S1.p4.1),[§2](https://arxiv.org/html/2607.00733#S2.p1.1),[§2](https://arxiv.org/html/2607.00733#S2.p2.1),[§3\.4](https://arxiv.org/html/2607.00733#S3.SS4.p1.1),[§4\.2](https://arxiv.org/html/2607.00733#S4.SS2.p1.3)\.
- R\. Tanese \(1989\)Distributed genetic algorithms for function optimization\.University of Michigan\.Cited by:[§3\.4](https://arxiv.org/html/2607.00733#S3.SS4.p1.1)\.
- L\. Wang, W\. Xu, Y\. Lan, Z\. Hu, Y\. Lan, R\. K\. Lee, and E\. Lim \(2023\)Plan\-and\-solve prompting: improving zero\-shot chain\-of\-thought reasoning by large language models\.InProceedings of the 61st annual meeting of the association for computational linguistics \(volume 1: long papers\),pp\. 2609–2634\.Cited by:[§2](https://arxiv.org/html/2607.00733#S2.p1.1)\.
- S\. N\. Wood \(2010\)Statistical inference for noisy nonlinear ecological dynamic systems\.Nature466\(7310\),pp\. 1102–1104\.Cited by:[§2](https://arxiv.org/html/2607.00733#S2.p3.1),[§3\.3](https://arxiv.org/html/2607.00733#S3.SS3.p2.12)\.
- J\. Yang, O\. Venkatachalam, M\. Kianezhad, S\. Vadgama, and R\. Yu \(2026\)Think like a scientist: physics\-guided llm agent for equation discovery\.arXiv preprint arXiv:2602\.12259\.Cited by:[§1](https://arxiv.org/html/2607.00733#S1.p3.1),[§2](https://arxiv.org/html/2607.00733#S2.p2.1)\.
- S\. Yao, D\. Yu, J\. Zhao, I\. Shafran, T\. Griffiths, Y\. Cao, and K\. Narasimhan \(2023\)Tree of thoughts: deliberate problem solving with large language models\.Advances in neural information processing systems36,pp\. 11809–11822\.Cited by:[§2](https://arxiv.org/html/2607.00733#S2.p1.1)\.
- S\. Yao, J\. Zhao, D\. Yu, N\. Du, I\. Shafran, K\. R\. Narasimhan, and Y\. Cao \(2022\)React: synergizing reasoning and acting in language models\.InThe eleventh international conference on learning representations,Cited by:[§1](https://arxiv.org/html/2607.00733#S1.p4.1),[§2](https://arxiv.org/html/2607.00733#S2.p1.1)\.

## Appendix

## Appendix ASummary Statistics

### A\.1Statistics for Within\-Structure Scoring \(logSL\)

We provide the full specification of the logSL summary statistic vector introduced in Section[3\.2](https://arxiv.org/html/2607.00733#S3.SS2)\. The logSL summary statistic vector is designed to be compact and interpretable, capturing the key distributional, temporal, and cross\-biomarker properties of the observed population while remaining computationally efficient\. It comprises the following four categories:

1. 1\.Marginal Distribution Statistics\.For each biomarker, all simulated values across all patients and time points are pooled and standardized using the mean and standard deviation of the observed data: z=y−μobsσobs\.z=\\frac\{y\-\\mu\_\{\\text\{obs\}\}\}\{\\sigma\_\{\\text\{obs\}\}\}\.Let\{z1,z2,…,zn\}\\\{z\_\{1\},z\_\{2\},\\ldots,z\_\{n\}\\\}denote the pooled standardized values\. After sorting, z\(1\)≤z\(2\)≤⋯≤z\(n\),qj=jn\+1,j=1,…,n,z\_\{\(1\)\}\\leq z\_\{\(2\)\}\\leq\\cdots\\leq z\_\{\(n\)\},\\qquad q\_\{j\}=\\frac\{j\}\{n\+1\},\\;j=1,\\dots,n,we approximate the empirical quantile function using a first\-degree polynomial: z\(j\)≈β0\+β1​qj\.z\_\{\(j\)\}\\approx\\beta\_\{0\}\+\\beta\_\{1\}q\_\{j\}\.The two coefficients\(β0,β1\)\(\\beta\_\{0\},\\beta\_\{1\}\)are included per biomarker, whereβ0\\beta\_\{0\}captures location shifts relative to the observed distribution andβ1\\beta\_\{1\}captures scale and tail spread\. This yields2​p2pstatistics in total\.
2. 2\.Temporal Dependence Statistics\.We compute lag\-1 Spearman autocorrelation for each biomarker\. Pooling all valid transitions across patients: ρ1=Corr​\(yt,yt−1\)\.\\rho\_\{1\}=\\mathrm\{Corr\}\(y\_\{t\},\\;y\_\{t\-1\}\)\.One autocorrelation coefficient is included per biomarker, yieldingppstatistics in total\.
3. 3\.Population\-level Trend Statistics\.For each biomarker, we compute the Spearman correlation between time and the population mean trajectory, capturing the overall directional trend of the biomarker over the observation window\. One trend statistic is included per biomarker, yieldingppstatistics in total\.
4. 4\.Dynamic Cross\-biomarker Statistics\.For biomarkersXXandYY, first differences are defined as: Δ​Xt=Xt−Xt−1,Δ​Yt=Yt−Yt−1\.\\Delta X\_\{t\}=X\_\{t\}\-X\_\{t\-1\},\\qquad\\Delta Y\_\{t\}=Y\_\{t\}\-Y\_\{t\-1\}\.Pooling all valid first differences across patients, we compute the Spearman correlation: ρΔ​X,Δ​Y=Corr​\(Δ​Xt,Δ​Yt\)\.\\rho\_\{\\Delta X,\\,\\Delta Y\}=\\mathrm\{Corr\}\(\\Delta X\_\{t\},\\;\\Delta Y\_\{t\}\)\.One correlation is included for each unordered biomarker pair, yieldingp​\(p−1\)2\\frac\{p\(p\-1\)\}\{2\}statistics in total\.

Forppbiomarkers, the resulting summary statistic vector has dimension4​p\+p​\(p−1\)24p\+\\frac\{p\(p\-1\)\}\{2\}\.

### A\.2Statistics for Across\-Structure Scoring and Diagnosis \(MNSD\)

We list the full set of statistics used for structure comparison via MNSD and as the diagnostic toolkit for the parameter inference agent \(Section[3\.5](https://arxiv.org/html/2607.00733#S3.SS5)\)\. Statistics with empirical IQR<0\.01<0\.01across observed trajectories are excluded from MNSD computation to avoid division by near\-zero values\.

Per\-variable statistics \(24 per biomarker\):

- •Distribution \(6\):mean, standard deviation, interquartile range, skewness, outlier fraction, value range\.
- •Autocorrelation \(6\):lag\-1, lag\-2, lag\-3 autocorrelation, first zero\-crossing of ACF, dominant frequency, spectral entropy\.
- •Stationarity \(2\):window mean variation, standard deviation of first differences\.
- •Entropy and Complexity \(2\):permutation entropy \(m=3m=3\), Lempel\-Ziv complexity\.
- •Scaling \(2\):detrended fluctuation analysis exponent, spectral slope\.
- •Model Fit \(3\):AR\(1\) coefficient, exponential smoothing alpha, turning point rate\.
- •Trend \(1\):Spearman correlation of values versus time\.
- •Volatility \(2\):mean crossing rate, mean absolute first difference\.

Cross\-variable statistics:

- •Level correlation:Spearman correlation of levels between two variables\.
- •Difference correlation:Spearman correlation of first differences\.
- •Cross\-correlation:cross\-correlation at lag 1\.

## Appendix BBaseline Implementation Details

We compare AgentODE against baselines encompassing a diverse range of methodologies, from sparse regression and evolutionary symbolic regression to LLM\-guided equation discovery and neural network\-based approaches\. Implementation details for each baseline are provided below\.

SINDySINDy\[Bruntonet al\.,[2016](https://arxiv.org/html/2607.00733#bib.bib16)\]is a sparse regression method that identifies governing equations by representing the time derivatives of state variables as a sparse linear combination of candidate library functions\. We use the PySINDy implementation\[De Silvaet al\.,[2020](https://arxiv.org/html/2607.00733#bib.bib45)\]222[https://github\.com/dynamicslab/pysindy](https://github.com/dynamicslab/pysindy)with a combined function library of polynomial terms up to degree 2 and nonlinear functions including exponential, logarithmic, trigonometric, and rational terms\. Derivatives are estimated using smoothed finite differences, and sparse regression is performed using the Sequentially Thresholded Least Squares \(STLSQ\) algorithm\. The sparsity threshold is tuned per system:λ=0\.02\\lambda=0\.02for Polymer DA and Apoptosis, andλ=0\.20\\lambda=0\.20for PKPD\-Immune\.

PySRPySR\[Cranmer,[2023](https://arxiv.org/html/2607.00733#bib.bib44)\]is an evolutionary symbolic regression method based on asynchronous multi\-island genetic programming\. We use the PySR package333[https://github\.com/MilesCranmer/PySR](https://github.com/MilesCranmer/PySR)and apply it by fitting one symbolic expression per state variable to model its time derivative, estimated from trajectories using Savitzky–Golay smoothing followed by finite differences\. The resulting expressions are combined into a full ODE system and integrated forward for evaluation\. We use 15 populations of size 33, with 1000 iterations per variable, and allow binary operators\{\+,−,×,÷\}\\\{\+,\-,\\times,\\div\\\}and unary operators\{exp,log,⋅,square\}\\\{\\exp,\\log,\\sqrt\{\\cdot\},\\text\{square\}\\\}with a maximum expression size of 20\.

LLM\-SR†We use the official implementation of LLM\-SR\[Shojaeeet al\.,[2024](https://arxiv.org/html/2607.00733#bib.bib5)\]444[https://github\.com/deep\-symbolic\-mathematics/LLM\-SR](https://github.com/deep-symbolic-mathematics/LLM-SR)and extend it from symbolic regression to ODE structure discovery\. Specifically, the LLM proposes candidate ODE structures, which are evaluated by simulating synthetic trajectories using numerical integration \(viaodeint\)\. For each proposed structure, parameters are fitted independently for each trajectory using BFGS to minimize trajectory\-level discrepancies\. This provides a reference setting with full access to individual\-level data\.

Neural ODENeural ODE\[Chenet al\.,[2018](https://arxiv.org/html/2607.00733#bib.bib32)\]learns system dynamics by parameterizing the right\-hand side of an ODE with a neural network,x˙=fθ​\(x\)\\dot\{x\}=f\_\{\\theta\}\(x\), and integrating it using a differentiable ODE solver\. We use thetorchdiffeqlibrary555[https://github\.com/rtqichen/torchdiffeq](https://github.com/rtqichen/torchdiffeq)with an MLP consisting of 3 hidden layers of size 64 andtanhactivations\. The model is trained using the Adam optimizer\[Kingma and Ba,[2014](https://arxiv.org/html/2607.00733#bib.bib46)\]with learning rate10−310^\{\-3\}and batch size 32 for up to 500 epochs with early stopping \(patience 50\) and learning rate decay\. Missing values are excluded from the loss computation\.

LLM\-based parameter initialization \(k=5k=5\)\.This ablation independently queries the LLMk=5k=5times with the ODE structure and problem context to propose parameter distributions, without visual summaries and without any diagnosis–update refinement\. The parameter distributions achieving the highest logSL score are selected\. This ablation serves as a lower bound, reflecting performance without any data\-driven iterative refinement\.

AgentODE w/o iterative refinement\.This ablation initializes parameter distributions using the parameter initialization prompt with visual summaries of the empirical data, identical to the full AgentODE framework\. However, the diagnosis–update loop is disabled and no iterative refinement is performed after initialization\. The resulting parameter distributions are used directly to score the candidate ODE structure via MNSD\.

## Appendix CAgentODE Implementation Details

LLM Backbone and ComputeAgentODE uses GPT\-5\.2 as the backbone LLM for both structure discovery and parameter inference, queried withmax\_tokens = 8192\. Temperature and top\-p are not explicitly set\. The primary computational cost comes from LLM API calls; ODE simulation supports both CPU and GPU execution and requires no specialized hardware\.

Structure DiscoveryThe outer loop maintains 10 islands\[Shojaeeet al\.,[2024](https://arxiv.org/html/2607.00733#bib.bib5)\], where each island independently evolves a population of candidate ODE structures to promote diversity, with periodic exchange of high\-performing structures across islands\. Each LLM call generates 4 candidate structures, with 2 high\-scoring structures fetched from the experience buffer to guide each proposal\. The buffer stores all structures grouped into clusters by−\-MNSD score\. Clusters are sampled via softmax over Z\-scored cluster scores, with temperature decaying from 0\.1 toward 0\.05 over 30,000 programs, favoring exploitation of high\-scoring structures while maintaining diversity\. Island resets occur every 4 hours, wiping the weakest half of islands and seeding each with the best structure from a surviving island\. A structure is registered in the buffer only if parameter inference yields at least one finite logSL, and is scored by−\-MNSD \(higher is better\) using a single simulation of 1,000 trajectories\.

logSL ComputationThe log synthetic likelihood is estimated usingNsim=150N\_\{\\text\{sim\}\}=150simulation replicates, each integrating 100 synthetic trajectories\. A trajectory is considered valid if all state variables are finite and fall within predefined physically plausible ranges\. A replicate is discarded if fewer than 80% of trajectories are valid, in which case logSL→−∞\\to\-\\infty\. A small regularization termλ=10−6\\lambda=10^\{\-6\}is added to the diagonal ofΣ^ψ\\hat\{\\Sigma\}\_\{\\psi\}before inversion to ensure numerical stability\.

Constraint violation handling\.When no parameter distributions proposed so far have yielded a finite logSL, the loop enters constraint violation mode\. A violation report is generated for the current parameter distributions, describing which state variables are out of physically plausible range and the fraction of invalid trajectories\. A dedicatedCONSTRAINT\_VIOLATIONprompt is constructed with the current parameters and violation report, and sent to the LLM to propose corrected distributions\.

Parameter InferenceThe inner diagnosis–update loop runs for a default of 10 iterations\. If the candidate structure yields a finite logSL and its−\-MNSD score exceeds the island’s current best, the loop is extended to 20 iterations, allowing more thorough parameter inference for promising structures without committing extra computation to poor ones\. Early stopping applies a patience of 5, counting both non\-improving steps and consecutive constraint\-violating steps\. The best iteration is selected by the highest logSL\. The diagnosis phase spans up to 5 tool\-calling turns, after which an update prompt is constructed with the 3 best parameter sets by logSL from the current iteration log\.

## Appendix DBenchmark Details

Standard benchmarks such as Lotka\-Volterra and SIR are well\-known to be susceptible to LLM recitation\. We therefore selected domain\-specific systems and verified their resistance to memorization through the contamination analysis in Table[3](https://arxiv.org/html/2607.00733#A4.T3)using the prompt provided in Appendix[G\.4](https://arxiv.org/html/2607.00733#A7.SS4)\.

Table 3:Contamination analysis for benchmark systems\. Recovery counts exact structural fingerprint matches across 5 runs usingGPT\-5\.2\.We provide the full ODE specifications, parameter settings, and data generation details for each benchmark below\.

### D\.1Apoptosis

The apoptosis system models programmed cell death through the interactions of three coupled protein concentrations: initiator protein, regulator protein, and effector protein\. The governing ODE system is adapted fromBideh and Gryak \[[2026](https://arxiv.org/html/2607.00733#bib.bib22)\]:

d​initiatord​t\\displaystyle\\frac\{d\\,\\text\{initiator\}\}\{dt\}=−k1⋅initiator⋅regulatorinitiator\+s1−k2⋅initiator\+k3\\displaystyle=\-\\frac\{k\_\{1\}\\cdot\\text\{initiator\}\\cdot\\text\{regulator\}\}\{\\text\{initiator\}\+s\_\{1\}\}\-k\_\{2\}\\cdot\\text\{initiator\}\+k\_\{3\}d​regulatord​t\\displaystyle\\frac\{d\\,\\text\{regulator\}\}\{dt\}=−k4⋅initiator⋅regulatorregulator\+s2−k5⋅regulatorregulator\+s3\+k6⋅effector⋅\(regulator\+s3\)\\displaystyle=\-\\frac\{k\_\{4\}\\cdot\\text\{initiator\}\\cdot\\text\{regulator\}\}\{\\text\{regulator\}\+s\_\{2\}\}\-\\frac\{k\_\{5\}\\cdot\\text\{regulator\}\}\{\\text\{regulator\}\+s\_\{3\}\}\+k\_\{6\}\\cdot\\text\{effector\}\\cdot\(\\text\{regulator\}\+s\_\{3\}\)d​effectord​t\\displaystyle\\frac\{d\\,\\text\{effector\}\}\{dt\}=k4⋅initiator⋅regulatorregulator\+s2\+k5⋅regulatorregulator\+s3−k6⋅effector⋅\(regulator\+s3\)\\displaystyle=\\phantom\{\-\}\\frac\{k\_\{4\}\\cdot\\text\{initiator\}\\cdot\\text\{regulator\}\}\{\\text\{regulator\}\+s\_\{2\}\}\+\\frac\{k\_\{5\}\\cdot\\text\{regulator\}\}\{\\text\{regulator\}\+s\_\{3\}\}\-k\_\{6\}\\cdot\\text\{effector\}\\cdot\(\\text\{regulator\}\+s\_\{3\}\)with fixed parametersk1=0\.40k\_\{1\}=0\.40,k2=0\.05k\_\{2\}=0\.05,k3=0\.10k\_\{3\}=0\.10,k4=7\.95k\_\{4\}=7\.95,k5=0\.20k\_\{5\}=0\.20,k6=0\.60k\_\{6\}=0\.60,s1=0\.10s\_\{1\}=0\.10,s2=2\.00s\_\{2\}=2\.00, ands3=0\.10s\_\{3\}=0\.10\. All parameters are fixed across individuals; only initial conditions vary\.

Data Generation\.A cohort of 100 trajectories is simulated overt∈\[0,10\]t\\in\[0,10\]with 101 evenly spaced time points \(Δ​t=0\.1\\Delta t=0\.1\)\. Initial conditions are sampled independently from uniform distributions: initiator∼𝒰​\(0\.001,0\.05\)\\sim\\mathcal\{U\}\(0\.001,0\.05\), regulator∼𝒰​\(0\.05,0\.50\)\\sim\\mathcal\{U\}\(0\.05,0\.50\), effector∼𝒰​\(0\.50,3\.00\)\\sim\\mathcal\{U\}\(0\.50,3\.00\)\. The cohort is split into 80 inference trajectories and 20 test trajectories drawn from the same initial condition distributions\.

### D\.2Polymer DA Cross\-linking

The Diels\-Alder \(DA\) cross\-linking system models the reversible covalent reaction between a diene and a dienophile component, a thermally reversible process widely studied in polymer chemistry for self\-healing and recyclable materials\. The system involves four coupled state variables: diene concentration \(x0x\_\{0\}\), dienophile concentration \(x1x\_\{1\}\), cross\-linked adduct concentration \(x2x\_\{2\}\), and de\-cross\-linked component concentration \(x3x\_\{3\}\)\. The governing ODE system is adapted fromBideh and Gryak \[[2026](https://arxiv.org/html/2607.00733#bib.bib22)\]:

d​x0d​t\\displaystyle\\frac\{d\\,x\_\{0\}\}\{dt\}=−k1​x0​x1−k2​x0​x3\+k3​x1\+k4​x2\\displaystyle=\-k\_\{1\}x\_\{0\}x\_\{1\}\-k\_\{2\}x\_\{0\}x\_\{3\}\+k\_\{3\}x\_\{1\}\+k\_\{4\}x\_\{2\}d​x1d​t\\displaystyle\\frac\{d\\,x\_\{1\}\}\{dt\}=−k1​x0​x1\+k2​x0​x3−k3​x1\+k4​x2\\displaystyle=\-k\_\{1\}x\_\{0\}x\_\{1\}\+k\_\{2\}x\_\{0\}x\_\{3\}\-k\_\{3\}x\_\{1\}\+k\_\{4\}x\_\{2\}d​x2d​t\\displaystyle\\frac\{d\\,x\_\{2\}\}\{dt\}=k1​x0​x1−k4​x2\\displaystyle=k\_\{1\}x\_\{0\}x\_\{1\}\-k\_\{4\}x\_\{2\}d​x3d​t\\displaystyle\\frac\{d\\,x\_\{3\}\}\{dt\}=−k2​x0​x3\+k3​x1\\displaystyle=\-k\_\{2\}x\_\{0\}x\_\{3\}\+k\_\{3\}x\_\{1\}
with fixed parametersk1=0\.03k\_\{1\}=0\.03,k2=0\.01k\_\{2\}=0\.01,k3=0\.08k\_\{3\}=0\.08, andk4=0\.03k\_\{4\}=0\.03\. All rate constants are fixed across individuals; only initial conditions vary\.

Data Generation\.The original system fromBideh and Gryak \[[2026](https://arxiv.org/html/2607.00733#bib.bib22)\]uses fixed initial conditions\. To simulate a heterogeneous population suitable for summary\-level inference, we extend the initial conditions to be sampled independently from uniform distributions:x0∼𝒰​\(1\.0,30\.0\)x\_\{0\}\\sim\\mathcal\{U\}\(1\.0,30\.0\),x1∼𝒰​\(0\.0,0\.1\)x\_\{1\}\\sim\\mathcal\{U\}\(0\.0,0\.1\),x2∼𝒰​\(0\.0,0\.1\)x\_\{2\}\\sim\\mathcal\{U\}\(0\.0,0\.1\),x3∼𝒰​\(0\.5,10\.0\)x\_\{3\}\\sim\\mathcal\{U\}\(0\.5,10\.0\), A cohort of 100 trajectories is simulated overt∈\[0,50\]t\\in\[0,50\]with 101 evenly spaced time points \(Δ​t=0\.5\\Delta t=0\.5\)\. The cohort is split into 80 inference trajectories and 20 test trajectories drawn from the same initial condition distributions\.

### D\.3PKPD\-Immune

The PKPD\-Immune system is a composite four\-variable ODE coupling two\-compartment pharmacokinetics with an immune suppression and effector cell response module\. The state variables are: drug amount in the absorption compartment \(depot\), drug amount in the central compartment \(central\), immune cell population \(immune\), and effector cell population \(effector\)\. The governing ODE system is:

d​depotd​t\\displaystyle\\frac\{d\\,\\text\{depot\}\}\{dt\}=−ka⋅depot\\displaystyle=\-k\_\{a\}\\cdot\\text\{depot\}d​centrald​t\\displaystyle\\frac\{d\\,\\text\{central\}\}\{dt\}=ka⋅depot−ke⋅central\\displaystyle=k\_\{a\}\\cdot\\text\{depot\}\-k\_\{e\}\\cdot\\text\{central\}d​immuned​t\\displaystyle\\frac\{d\\,\\text\{immune\}\}\{dt\}=s−d⋅immune−k1​\(central/V\)⋅immuneI​C50\+central/V\\displaystyle=s\-d\\cdot\\text\{immune\}\-\\frac\{k\_\{1\}\(\\text\{central\}/V\)\\cdot\\text\{immune\}\}\{IC\_\{50\}\+\\text\{central\}/V\}d​effectord​t\\displaystyle\\frac\{d\\,\\text\{effector\}\}\{dt\}=ρ⋅immune⋅effectorimmune\+σ−δ⋅effector−ε⋅effector2\\displaystyle=\\frac\{\\rho\\cdot\\text\{immune\}\\cdot\\text\{effector\}\}\{\\text\{immune\}\+\\sigma\}\-\\delta\\cdot\\text\{effector\}\-\\varepsilon\\cdot\\text\{effector\}^\{2\}whereε=0\.05\\varepsilon=0\.05is a fixed structural parameter governing effector intraspecific competition\. Both initial conditions and model parameters vary across individuals, reflecting pharmacokinetic and pharmacodynamic variability in drug response\. All parameters are sampled independently from log\-normal distributions; Table[4](https://arxiv.org/html/2607.00733#A4.T4)reports the mean and standard deviation in linear space\.

Table 4:PKPD\-Immune parameter distributions \(log\-normal\)\.Data Generation\.A cohort of 100 trajectories is simulated overt∈\[0,50\]t\\in\[0,50\]days with 101 evenly spaced time points \(Δ​t=0\.5\\Delta t=0\.5days\)\. Initial conditions are sampled from log\-normal distributions: depot∼LogNormal​\(μ=5\.0,σ=1\.0\)\\sim\\text\{LogNormal\}\(\\mu=5\.0,\\sigma=1\.0\), central=0=0\(central compartment always empty at dose time\), immune∼LogNormal​\(μ=6\.67,σ=1\.5\)\\sim\\text\{LogNormal\}\(\\mu=6\.67,\\sigma=1\.5\)initialized near the drug\-free immune steady state, and effector∼LogNormal​\(μ=1\.0,σ=0\.30\)\\sim\\text\{LogNormal\}\(\\mu=1\.0,\\sigma=0\.30\)\. The cohort is split into 80 inference trajectories and 20 test trajectories drawn from the same distributions\.

## Appendix EClinical Dataset Details

### E\.1Acute Kidney Injury \(AKI\)

Patient Selection\.Patients aged 50–64 years at admission are selected from MIMIC\-IV 3\.1\[Johnsonet al\.,[2023](https://arxiv.org/html/2607.00733#bib.bib39)\]666[https://physionet\.org/content/mimiciv/3\.1/](https://physionet.org/content/mimiciv/3.1/)\. AKI is confirmed using ICD\-10 diagnosis codes \(N17\.\*\) and/or KDIGO clinical criteria \(≥\\geq0\.3 mg/dL creatinine increase within 48 hours, or≥\\geq1\.5×\\timesbaseline within 7 days\)\. Only admissions with at least 20 unique timepoints within the first 168 hours and at least two of the three core biomarkers present at each timepoint are retained\. The final cohort comprises 345 patients \(353 admissions, 8,249 timepoints; mean age57\.4±4\.457\.4\\pm 4\.4years\)\.

Train/Test Split\.The cohort is split at the admission level into 80% inference \(282 patients\) and 20% test \(71 patients\)\.

Preprocessing\.No interpolation or imputation is applied\. Missing values and irregular sampling are handled natively by AgentODE through population\-level summary statistics, and by Neural ODE and LLM\-SR\-ODE through masked loss computation\.

Ethics\.MIMIC\-IV is accessed under the PhysioNet Credentialed Health Data License 1\.5\.0 and the PhysioNet Credentialed Health Data Use Agreement 1\.5\.0\. No individual\-level data are released as part of this work\.

Preprocessing\.No interpolation or imputation is applied to the raw measurements\. Missing values and irregular sampling are handled natively by AgentODE through population\-level summary statistics, and by Neural ODE and LLM\-SR\-ODE through masked loss computation\.

### E\.2Recessive Dystrophic Epidermolysis Bullosa \(RDEB\)

Dataset\.The dataset comprises 46 patients with severe RDEB\[Bardhanet al\.,[2020](https://arxiv.org/html/2607.00733#bib.bib42), Reimeret al\.,[2020](https://arxiv.org/html/2607.00733#bib.bib40), Karakioulakiet al\.,[2026](https://arxiv.org/html/2607.00733#bib.bib41)\]\. Patients ranged from 0 to 30\.4 years of age at their first recorded visit\. We focus on longitudinal follow\-up spanningt∈\[0,144\]t\\in\[0,144\]months since first hospital visit \(mean 53\.3 months per patient\)\. A total of 231 observations are available across 46 patients, with a mean of 5\.0 timepoints per patient and 16 patients having as few as 2 timepoints\.

State Variables\.Four state variables are modeled: serum albumin \(g/dL\), log\-transformed C\-reactive protein \(log10\(CRP\+\+1\), where CRP is in mg/L\), hemoglobin \(g/dL\), and BMI relative to age\- and sex\-matched WHO reference \(BMIrel, kg/m2\)\. BMIrelis computed as the deviation of each patient’s BMI from the age\- and sex\-matched WHO healthy reference, such that a value of 0 indicates a healthy BMI and negative values indicate underweight relative to peers\.

Preprocessing\.CRP is log10\(CRP\+\+1\) transformed to stabilize variance\. BMIrelis computed as the deviation from age\- and sex\-matched WHO reference values\. No interpolation or imputation is applied\.

Evaluation\.Due to the limited cohort size, no train/test split is performed and results are reported in\-sample, as noted in Table[2](https://arxiv.org/html/2607.00733#S4.T2)\.

Ethics\.The RDEB dataset was collected under institutional ethics approval at the authors’ institution\. All patients or their legal guardians provided informed consent\. No individual\-level data are released as part of this work\.

## Appendix FAdditional Results

### F\.1Apoptosis

Figure[4](https://arxiv.org/html/2607.00733#A6.F4)compares the ODE structures recovered by different methods on the Apoptosis benchmark\. SINDy recovers highly complex equations containing many redundant nonlinear, trigonometric, and exponential terms, producing expressions that are difficult to interpret mechanistically despite achieving strong numerical performance\. PySR discovers substantially more compact equations, but the recovered interactions deviate from the underlying system structure and fail to preserve the coupled regulatory dynamics between variables\. LLM\-SR†recovers a sparse and biologically interpretable structure that closely resembles the ground\-truth interaction pattern, correctly identifying key cross\-variable couplings between initiator, regulator, and effector dynamics\. AgentODE similarly recovers a compact interaction structure\. Although the exact symbolic form differs from the ground truth, the discovered system preserves the core dynamical relationships between variables, including inhibitory and recruitment\-like interactions that govern the population\-level apoptosis behavior\.

![Refer to caption](https://arxiv.org/html/2607.00733v1/x4.png)Figure 4:Recovered ODE structures for the Apoptosis benchmark using AgentODE, LLM\-SR†, SINDy, and PySR\.
### F\.2Polymer DA Cross\-linking

Figure[5](https://arxiv.org/html/2607.00733#A6.F5)compares the ODE structures recovered by different methods on the Polymer DA Cross\-linking benchmark\. The ground\-truth system is governed by sparse reaction\-like terms involving cross\-linking between reactants, product formation, and de\-cross\-linking dynamics\. SINDy again recovers highly complex equations with many nonlinear, exponential, logarithmic, and trigonometric terms, making the resulting system difficult to interpret mechanistically\. PySR produces more compact equations, but the recovered expressions contain nonphysical functional forms such as logarithmic and high\-order rational terms, which do not align well with the underlying reaction kinetics\. LLM\-SR†recovers an interpretable reversible reaction network, but its proposed structure differs from the ground\-truth reaction topology\. In contrast, AgentODE recovers a compact mass\-action\-like system with forward cross\-linking, reverse de\-cross\-linking, and side\-interaction terms\. Although the exact coefficients and terms differ from the ground truth, the recovered structure preserves the main mechanistic motifs of the Polymer system, including reactant\-product conversion and coupling through the de\-cross\-linked component\.

![Refer to caption](https://arxiv.org/html/2607.00733v1/x5.png)Figure 5:Recovered ODE structures for the Polymer DA Cross\-linking benchmark using AgentODE, LLM\-SR†, SINDy, and PySR\.
### F\.3PKPD\-Immune

Figure[6](https://arxiv.org/html/2607.00733#A6.F6)compares the ODE structures recovered by different methods on the PKPD\-Immune benchmark\. The ground\-truth system combines two\-compartment pharmacokinetics with immune suppression and effector\-cell recruitment dynamics\. SINDy recovers a simplified structure that captures parts of the depot and central compartment dynamics, but fails to model effector dynamics and omits key immune–effector interactions\. PySR produces compact equations but introduces nonphysical terms, including logarithmic effector feedback, and does not recover the mechanistic immune recruitment structure\. In contrast, both LLM\-SR†and AgentODE recover interpretable PKPD\-style systems\. LLM\-SR†closely matches the ground\-truth structure, including absorption, elimination, saturable immune suppression, and effector decay\. AgentODE recovers the same core pharmacokinetic and immune\-suppression mechanisms while using only summary\-level statistics\. Although its effector recruitment term differs from the ground truth, the recovered structure preserves the main causal pathway from depot absorption to central drug exposure, immune suppression, and effector dynamics\.

![Refer to caption](https://arxiv.org/html/2607.00733v1/x6.png)Figure 6:Recovered ODE structures for the PKPD\-Immune benchmark using AgentODE, LLM\-SR†, SINDy, and PySR\.
### F\.4AKI

#### ODE structure discovered by AgentODE\.

AgentODE discovers a linear coupled biomarker system:

d​Crd​t\\displaystyle\\frac\{d\\mathrm\{Cr\}\}\{dt\}=θ0−θ1​Cr\+θ2​BUN,\\displaystyle=\\theta\_\{0\}\-\\theta\_\{1\}\\mathrm\{Cr\}\+\\theta\_\{2\}\\mathrm\{BUN\},d​BUNd​t\\displaystyle\\frac\{d\\mathrm\{BUN\}\}\{dt\}=θ3−θ4​BUN\+θ5​Cr,\\displaystyle=\\theta\_\{3\}\-\\theta\_\{4\}\\mathrm\{BUN\}\+\\theta\_\{5\}\\mathrm\{Cr\},d​Kd​t\\displaystyle\\frac\{d\\mathrm\{K\}\}\{dt\}=θ6−θ7​K\+θ8​Cr\+θ9​BUN\.\\displaystyle=\\theta\_\{6\}\-\\theta\_\{7\}\\mathrm\{K\}\+\\theta\_\{8\}\\mathrm\{Cr\}\+\\theta\_\{9\}\\mathrm\{BUN\}\.
This structure represents each biomarker through a baseline source term, a clearance or recovery term, and cross\-biomarker coupling\. Creatinine and BUN are mutually coupled, reflecting shared renal dysfunction and azotemia dynamics, while potassium is driven by both creatinine and BUN together with its own elimination term\. The resulting system is parsimonious and clinically interpretable\.

#### ODE structure discovered by LLM\-SR†\.

LLM\-SR†discovers the following AKI system:

d​Crd​t\\displaystyle\\frac\{d\\mathrm\{Cr\}\}\{dt\}=θ0\+θ1​Cr−θ2​Cr,\\displaystyle=\\theta\_\{0\}\+\\theta\_\{1\}\\mathrm\{Cr\}\-\\theta\_\{2\}\\mathrm\{Cr\},d​BUNd​t\\displaystyle\\frac\{d\\mathrm\{BUN\}\}\{dt\}=θ3\+θ4​Cr\+θ5​BUN−θ6​BUN,\\displaystyle=\\theta\_\{3\}\+\\theta\_\{4\}\\mathrm\{Cr\}\+\\theta\_\{5\}\\mathrm\{BUN\}\-\\theta\_\{6\}\\mathrm\{BUN\},d​Kd​t\\displaystyle\\frac\{d\\mathrm\{K\}\}\{dt\}=θ7​Cr\+θ8​BUN\+θ9−θ10​K\.\\displaystyle=\\theta\_\{7\}\\mathrm\{Cr\}\+\\theta\_\{8\}\\mathrm\{BUN\}\+\\theta\_\{9\}\-\\theta\_\{10\}\\mathrm\{K\}\.
This structure also captures clinically plausible AKI dynamics, including creatinine\-driven BUN changes and potassium dependence on renal dysfunction markers\. However, the positive and negative self\-effects for creatinine and BUN reduce to net linear turnover terms, making these effects harder to distinguish under sparse clinical sampling\. Compared with AgentODE, the LLM\-SR†structure is slightly less parsimonious and lacks direct BUN\-to\-creatinine coupling\.

### F\.5RDEB

#### ODE structure developed independently by ODE experts in collaboration with RDEB clinicians\.

The expert\-derived model uses four state variables: nutrition \(NN\), inflammation \(II\), anemia state \(AA\), and general health \(HH\), along with an unobserved auxiliary inflammatory driveI~\\widetilde\{I\}\. These states were mapped to the previously used clinical proxy variables albumin, C\-reactive protein, hemoglobin, and BMI relative to age\- and sex\-matched WHO reference values, respectively\. To facilitate quantitative interpretation of model parameters, clinical variables were standardized to a\[0,1\]\[0,1\]range\. The model is specified by:

d​Nd​t\\displaystyle\\frac\{dN\}\{dt\}=λN​\(N∗​\(I\)−N\),N∗​\(I\)=kI,N​\(1−I\)kI,N\+I\\displaystyle=\\lambda\_\{N\}\\left\(N^\{\*\}\(I\)\-N\\right\),\\quad N^\{\*\}\(I\)=\\frac\{k\_\{I,N\}\(1\-I\)\}\{k\_\{I,N\}\+I\}d​Id​t\\displaystyle\\frac\{dI\}\{dt\}=λI​\(I∗​\(H,I~\)−I\),I∗​\(H,I~\)=11\+I~​\[kH,I​\(1−H\)kH,I\+H\+I~\]\\displaystyle=\\lambda\_\{I\}\\left\(I^\{\*\}\(H,\\widetilde\{I\}\)\-I\\right\),\\quad I^\{\*\}\(H,\\widetilde\{I\}\)=\\frac\{1\}\{1\+\\widetilde\{I\}\}\\left\[\\frac\{k\_\{H,I\}\(1\-H\)\}\{k\_\{H,I\}\+H\}\+\\widetilde\{I\}\\right\]d​Ad​t\\displaystyle\\frac\{dA\}\{dt\}=λA​\(A∗​\(I\)−A\),A∗​\(I\)=kI,A​\(1−I\)kI,A\+I\\displaystyle=\\lambda\_\{A\}\\left\(A^\{\*\}\(I\)\-A\\right\),\\quad A^\{\*\}\(I\)=\\frac\{k\_\{I,A\}\(1\-I\)\}\{k\_\{I,A\}\+I\}d​Hd​t\\displaystyle\\frac\{dH\}\{dt\}=λH​\(H∗​\(N,A,I\)−H\),H∗​\(N,A,I\)=13​\(kI,H​\(1−I\)kI,H\+I\+\(1\+kN,H\)​NkN,H\+N\+\(1\+kA,H\)​AkA,H\+A\)\\displaystyle=\\lambda\_\{H\}\\left\(H^\{\*\}\(N,A,I\)\-H\\right\),\\quad H^\{\*\}\(N,A,I\)=\\frac\{1\}\{3\}\\left\(\\frac\{k\_\{I,H\}\(1\-I\)\}\{k\_\{I,H\}\+I\}\+\\frac\{\(1\+k\_\{N,H\}\)N\}\{k\_\{N,H\}\+N\}\+\\frac\{\(1\+k\_\{A,H\}\)A\}\{k\_\{A,H\}\+A\}\\right\)d​I~d​t\\displaystyle\\frac\{d\\widetilde\{I\}\}\{dt\}=pI,I~​\(0\)=0\\displaystyle=p\_\{I\},\\quad\\widetilde\{I\}\(0\)=0
where each observed state relaxes toward a state\-specific set point that depends on the current value of the other states\. The auxiliary stateI~\\widetilde\{I\}represents the unobserved progressive inflammatory drive and is not used as an observed clinical variable\.

#### ODE structure discovered by AgentODE\.

The following ODE structure models the coupled dynamics of serum albumin, log\-transformed C\-reactive protein \(log\\log\-CRP\), hemoglobin, and BMI relative to age\- and sex\-matched WHO reference \(BMIrel\)\.

d​albumind​t\\displaystyle\\frac\{d\\,\\text\{albumin\}\}\{dt\}=θ0​\(θ1−albumin\)−θ2⋅log⁡CRP\\displaystyle=\\theta\_\{0\}\\left\(\\theta\_\{1\}\-\\text\{albumin\}\\right\)\-\\theta\_\{2\}\\cdot\\log\\text\{CRP\}d​log⁡CRPd​t\\displaystyle\\frac\{d\\,\\log\\text\{CRP\}\}\{dt\}=θ3​\(θ4−log⁡CRP\)\+θ5​\(θ1−albumin\)\\displaystyle=\\theta\_\{3\}\\left\(\\theta\_\{4\}\-\\log\\text\{CRP\}\\right\)\+\\theta\_\{5\}\\left\(\\theta\_\{1\}\-\\text\{albumin\}\\right\)d​Hbd​t\\displaystyle\\frac\{d\\,\\text\{Hb\}\}\{dt\}=θ6​\(θ7−Hb\)−θ8⋅log⁡CRP\+θ9​\(albumin−θ1\)\\displaystyle=\\theta\_\{6\}\\left\(\\theta\_\{7\}\-\\text\{Hb\}\\right\)\-\\theta\_\{8\}\\cdot\\log\\text\{CRP\}\+\\theta\_\{9\}\\left\(\\text\{albumin\}\-\\theta\_\{1\}\\right\)d​BMIreld​t\\displaystyle\\frac\{d\\,\\text\{BMI\}\_\{\\text\{rel\}\}\}\{dt\}=θ10​\(θ11−BMIrel\)\+θ92​\(albumin−θ1\)\\displaystyle=\\theta\_\{10\}\\left\(\\theta\_\{11\}\-\\text\{BMI\}\_\{\\text\{rel\}\}\\right\)\+\\frac\{\\theta\_\{9\}\}\{2\}\\left\(\\text\{albumin\}\-\\theta\_\{1\}\\right\)
Comparison\.The expert\-derived model and the model found by AgentODE approach the phenomenological relationships between variables differently, but both recognize the key disease loop, as visualized in Figure[7](https://arxiv.org/html/2607.00733#A6.F7): inflammation \(log\-CRP\) drives worsening of the other clinical variables, while a second variable attenuates inflammation\. In the expert model, BMIrel\{\}\_\{\\text\{rel\}\}plays this protective role as a proxy for patient health, while AgentODE identifies albumin as the protective variable, interpreted as a proxy for nutritional status\. Both models plausibly link the clinical variables into a phenomenological description, for which no unique solution exists\.

#### ODE structure discovered by AgentODE w/o iterative refinement\.

The following structure achieves a marginally lower RMSE \(0\.628 vs 0\.644\) than the full AgentODE model\. However, mechanistic inspection reveals several implausibilities\. The BMIrel\{\}\_\{\\text\{rel\}\}equation reuses parametersθ0,θ1,θ2,θ3\\theta\_\{0\},\\theta\_\{1\},\\theta\_\{2\},\\theta\_\{3\}from the albumin equation, imposing the constraint that both variables share identical relaxation rates and baselines, which is clinically implausible as they reflect distinct physiological processes\. More fundamentally, higher BMIrel\{\}\_\{\\text\{rel\}\}negatively affects albumin and hemoglobin, whereas higher BMIrel\{\}\_\{\\text\{rel\}\}generally reflects better nutritional status and should be associated with improved clinical markers\. Conversely, higher albumin negatively affects BMIrel\{\}\_\{\\text\{rel\}\}, which contradicts the established association between good nutritional status and healthy body weight\. As shown in Figure[7](https://arxiv.org/html/2607.00733#A6.F7), the resulting causal structure is dense with exclusively negative cross\-variable effects and no clear protective mechanism, in contrast to the parsimonious and clinically interpretable structure recovered by AgentODE\. This illustrates that in sparse and noisy clinical settings, lower RMSE may reflect overfitting to noise rather than a more mechanistically valid model\.

d​albumind​t\\displaystyle\\frac\{d\\,\\text\{albumin\}\}\{dt\}=θ0​\(θ1−albumin\)−θ2⋅log⁡CRP−θ3⋅BMIrel\\displaystyle=\\theta\_\{0\}\\left\(\\theta\_\{1\}\-\\text\{albumin\}\\right\)\-\\theta\_\{2\}\\cdot\\log\\text\{CRP\}\-\\theta\_\{3\}\\cdot\\text\{BMI\}\_\{\\text\{rel\}\}d​log⁡CRPd​t\\displaystyle\\frac\{d\\,\\log\\text\{CRP\}\}\{dt\}=θ4​\(θ5−log⁡CRP\)−θ6⋅albumin−θ7⋅BMIrel\\displaystyle=\\theta\_\{4\}\\left\(\\theta\_\{5\}\-\\log\\text\{CRP\}\\right\)\-\\theta\_\{6\}\\cdot\\text\{albumin\}\-\\theta\_\{7\}\\cdot\\text\{BMI\}\_\{\\text\{rel\}\}d​Hbd​t\\displaystyle\\frac\{d\\,\\text\{Hb\}\}\{dt\}=θ8​\(θ9−Hb\)−θ10⋅log⁡CRP−θ11⋅BMIrel\\displaystyle=\\theta\_\{8\}\\left\(\\theta\_\{9\}\-\\text\{Hb\}\\right\)\-\\theta\_\{10\}\\cdot\\log\\text\{CRP\}\-\\theta\_\{11\}\\cdot\\text\{BMI\}\_\{\\text\{rel\}\}d​BMIreld​t\\displaystyle\\frac\{d\\,\\text\{BMI\}\_\{\\text\{rel\}\}\}\{dt\}=θ0​\(θ1−BMIrel\)−θ2⋅log⁡CRP−θ3⋅albumin\\displaystyle=\\theta\_\{0\}\\left\(\\theta\_\{1\}\-\\text\{BMI\}\_\{\\text\{rel\}\}\\right\)\-\\theta\_\{2\}\\cdot\\log\\text\{CRP\}\-\\theta\_\{3\}\\cdot\\text\{albumin\}

#### ODE structure discovered by LLM\-SR†\.

This structure achieves a lower RMSE \(0\.532\) than AgentODE \(0\.644\) by treating BMIrel\{\}\_\{\\text\{rel\}\}as a static variable with zero derivative\. While this trivially reduces prediction error on a sparse and noisy variable, it entirely forgoes modeling BMIrel\{\}\_\{\\text\{rel\}\}dynamics\. Since BMIrel\{\}\_\{\\text\{rel\}\}is a key indicator of disease progression and nutritional status in RDEB, a model that does not capture its temporal evolution cannot provide mechanistic insights into this clinically important dimension\. Furthermore, the albumin and log\-CRP equations use constant baseline synthesis terms \(θ0\\theta\_\{0\}andθ4\\theta\_\{4\}\) rather than relaxation toward a physiologically meaningful baseline, which limits biological interpretability\. As shown in Figure[7](https://arxiv.org/html/2607.00733#A6.F7), BMIrel\{\}\_\{\\text\{rel\}\}is entirely disconnected from the causal structure\. This further illustrates that lower RMSE does not necessarily reflect a more mechanistically informative model in sparse and noisy clinical settings\.

d​albumind​t\\displaystyle\\frac\{d\\,\\text\{albumin\}\}\{dt\}=θ0−θ1⋅albumin−θ2⋅log⁡CRP−θ3⋅BMIrel\\displaystyle=\\theta\_\{0\}\-\\theta\_\{1\}\\cdot\\text\{albumin\}\-\\theta\_\{2\}\\cdot\\log\\text\{CRP\}\-\\theta\_\{3\}\\cdot\\text\{BMI\}\_\{\\text\{rel\}\}d​log⁡CRPd​t\\displaystyle\\frac\{d\\,\\log\\text\{CRP\}\}\{dt\}=θ4−θ5⋅log⁡CRP−θ6⋅albumin−θ7⋅BMIrel\\displaystyle=\\theta\_\{4\}\-\\theta\_\{5\}\\cdot\\log\\text\{CRP\}\-\\theta\_\{6\}\\cdot\\text\{albumin\}\-\\theta\_\{7\}\\cdot\\text\{BMI\}\_\{\\text\{rel\}\}d​Hbd​t\\displaystyle\\frac\{d\\,\\text\{Hb\}\}\{dt\}=θ8−θ9⋅Hb−θ10⋅log⁡CRP\+θ11⋅albumin\\displaystyle=\\theta\_\{8\}\-\\theta\_\{9\}\\cdot\\text\{Hb\}\-\\theta\_\{10\}\\cdot\\log\\text\{CRP\}\+\\theta\_\{11\}\\cdot\\text\{albumin\}d​BMIreld​t\\displaystyle\\frac\{d\\,\\text\{BMI\}\_\{\\text\{rel\}\}\}\{dt\}=0\\displaystyle=0

#### Analysis\.

These comparisons highlight a fundamental advantage of summary\-level inference in clinical settings\. Methods that optimize directly against individual trajectories, such as LLM\-SR†, can achieve low RMSE by overfitting to noise rather than recovering the true underlying mechanism\. In contrast, reasoning from population\-level summary statistics encourages the LLM to identify patterns that are robust across the cohort, promoting mechanistically grounded structure discovery\. This distinction becomes particularly consequential in the experience buffer: since high\-scoring structures are retrieved to guide future proposals, a mechanistically flawed structure with artificially low RMSE can mislead the search toward implausible regions of the structure space\. Summary statistics, by averaging out individual\-level noise, provide a more representative and informative signal for clinical ODE discovery than trajectory\-level fitting alone\.

![Refer to caption](https://arxiv.org/html/2607.00733v1/x7.png)Figure 7:Causal diagrams for the RDEB ODE structures from \(first\) the expert\-derived, \(second\) AgentODE, \(third\) AgentODE w/o iterative refinement, and \(fourth\) LLM\-SR†\. Red arrows indicate negative effects, blue arrows indicate positive effects\.
#### Inferred parameter distributions

Figure[8](https://arxiv.org/html/2607.00733#A6.F8)compares the parameter distributions inferred by AgentODE against empirical distributions for the RDEB dataset\. To obtain the empirical distributions, we fit the discovered ODE structure to each individual trajectory independently via BFGS optimization, yielding one parameter vector per patient\. The resulting per\-patient parameter values are shown as histograms, with a log\-normal curve fitted to each, and compared against the log\-normal distributions inferred by AgentODE from population\-level summary statistics alone\.

For the baseline and relaxation parameters \(θ1,θ4,θ7,θ10\\theta\_\{1\},\\theta\_\{4\},\\theta\_\{7\},\\theta\_\{10\}\), the inferred and fitted distributions are in close agreement\. For the cross\-variable coupling parameters \(θ2,θ5,θ8,θ9\\theta\_\{2\},\\theta\_\{5\},\\theta\_\{8\},\\theta\_\{9\}\), the inferred and fitted distributions show certain discrepancies in magnitude, reflecting the inherent difficulty of precisely estimating interaction strengths in this setting: the data are sparse, heterogeneous, and subject to largely missing pairwise observations, only four clinical variables are modeled, and the mechanistic coupling between these variables in RDEB remains an open clinical question with limited established knowledge\. Nevertheless, Figure[9](https://arxiv.org/html/2607.00733#A6.F9)shows that the signs of most pairwise first\-difference Spearman correlations are correctly recovered, suggesting that AgentODE captures the directional structure of cross\-variable interactions even when magnitudes remain uncertain\. The exception is albumin–BMIrel\{\}\_\{\\text\{rel\}\}, which is empirically near zero \(−0\.006\-0\.006\) and carries insufficient signal\.

![Refer to caption](https://arxiv.org/html/2607.00733v1/x8.png)Figure 8:Parameter distribution comparisons for the AgentODE\-discovered RDEB ODE structure\.Blue histograms: per\-patient parameter values from BFGS fitting to individual trajectories\.Blue solid curves: log\-normal fits to these values\.Orange dashed curves: log\-normal distributions inferred by AgentODE from population\-level summary statistics\.![Refer to caption](https://arxiv.org/html/2607.00733v1/x9.png)Figure 9:RDEB first\-difference Spearman correlation heatmaps for empirical and synthetic trajectories from AgentODE discovered ODE structure\.

## Appendix GPrompts

All prompts in this section use the PKPD\-Immune benchmark as a representative example\. Variable names, units, time horizon, and domain context are substituted accordingly for other datasets\.

### G\.1Structure Discovery Prompt

Structure Discovery Prompt \(PKPD\-Immune\)Youareahelpfulassistanttaskedwithdiscoveringordinary

differentialequations\(ODEs\)structuresfordynamicprocesses\.

Completethe’system’functionbelow,consideringtherolesof

eachinput,theirinteractions,andplausibledynamics\.Let’s

thinkstepbystep,butprovideonly\(1\)afinalbriefexplanation

\(2\-5sentences\)ofyourreasoning,limitationsandunmodeled

effects,followedby\(2\)thecompletedODEexpressions\.

"""

FindthesymbolicskeletonofanODEsystemrepresenting

afour\-variablepharmacokinetic\-immunodynamicprocess\.

Thesysteminvolvesfourcoupledstatevariables:

\-depot\(amount\):Drugamountintheabsorptioncompartment\.

\-central\(amount\):Drugamountinthecentralcompartment\.

\-immune\(cells\):Immunecellpopulationsubjectto

drug\-mediatedsuppressionfollowingsaturatingdose\-response\.

\-effector\(cells\):Effectorcellpopulationrecruitedthrough

contact\-dependentinteractionwithimmunecells\.

Thetimeunitisdays\.Observationsspantin\[0,50\.0\]

HARDCONSTRAINT:

\-TheparameterarrayhasatmostMAX\_NPARAMS=10slots\.

\-DoNOTskipindices\.

\-Useonlyasmanyparametersasnecessary;avoidintroducing

redundantparameters\.

\-Youmustnotuseanyparams\[k\]withk\>=MAX\_NPARAMS\.

"""

importnumpyasnp

MAX\_NPARAMS=10

params=\[1\.0\]\*MAX\_NPARAMS

defsystem\_v0\(biomarkers:np\.ndarray,params:np\.ndarray\)\-\>np\.ndarray:

"""Symbolicskeletonfora4\-variableODEsystem\.

Args:

biomarkers:Anumpyarraycontaining

\[depot,central,immune,effector\]\.

params:Arrayofnumericconstantsorparameterstobeoptimized\.

Youmustnotuseanyparams\[k\]withk\>=MAX\_NPARAMS\.

Returns:

Anumpyarrayrepresentingthetemporalderivatives\(d/dt\)

ofthestatevariables

\[d\_depot/dt,d\_central/dt,d\_immune/dt,d\_effector/dt\]\.

"""

depot,central,immune,effector=biomarkers

\#\-\-\-\-depotdynamics\-\-\-\-

d\_depot=params\[0\]\*depot

\#\-\-\-\-centraldynamics\-\-\-\-

d\_central=params\[1\]\*central

\#\-\-\-\-immunedynamics\-\-\-\-

d\_immune=params\[2\]\*immune

\#\-\-\-\-effectordynamics\-\-\-\-

d\_effector=params\[3\]\*effector

d\_biomarkers=np\.array\(\[d\_depot,d\_central,d\_immune,d\_effector\],

dtype=float\)

returnd\_biomarkers

defsystem\_v1\(biomarkers:np\.ndarray,params:np\.ndarray\)\-\>np\.ndarray:

"""Improvedversionof‘system\_v0‘\."""

### G\.2Parameter Initialization Prompt

Youareanexpertdynamicalsystemsanalystwithdeepexpertisein

mechanisticmodelingandparameterdistributioninferencefor

coupledpharmacologicalandimmunologicalsystems\.

Yourtaskistoproposeandjustifymoderatelyinformativepriors

\(linear\-spacemeanandSD\)fortheparametersofthemechanistic

ODEmodelbelow,suchthatsimulatedtrajectoriesalignwiththe

observeddataovertin\[0,50\]\.

DEFINITION:MODERATELYINFORMATIVE

Apriorismoderatelyinformativeifit:

\-Rulesoutphysicallyimplausiblevalues

\-Isconsistentwiththeobservedtrajectorypatterns

\-Reflectsthedegreeofvariabilitypresentinthedata\-\-\-

neithertoodiffusenortooconcentrated

Simulation:tindays,horizon=50

Statevariables:depot\(drugdepotamount\),central\(drugcentral

compartmentamount\),immune\(immunecellpopulation\),effector

\(effectorcellpopulation\)

defsystem\(biomarkers:np\.ndarray,params:np\.ndarray\)\-\>np\.ndarray:

depot,central,immune,effector=biomarkers

\#params\[0\]=absorptionratefromdepottocentral

\#params\[1\]=centraleliminationrate

\#params\[2\]=drugeffectscaleonimmunesuppression

\#params\[3\]=drughalf\-saturationconstant

\#params\[4\]=immunebaseline/sourcerate

\#params\[5\]=immunenaturaldecayrate

\#params\[6\]=effectorrecruitmentratefromimmunecontact

\#params\[7\]=effectornaturaldecayrate

k\_a=params\[0\]

k\_el=params\[1\]

e\_max=params\[2\]

ec50=params\[3\]

s\_imm=params\[4\]

d\_imm=params\[5\]

k\_recruit=params\[6\]

d\_eff=params\[7\]

drug\_effect=e\_max\*central/\(ec50\+central\+1e\-12\)

d\_depot=\-k\_a\*depot

d\_central=k\_a\*depot\-k\_el\*central

d\_immune=s\_imm\-d\_imm\*immune\-drug\_effect\*immune

d\_effector=k\_recruit\*immune\*effector\-d\_eff\*effector

returnnp\.array\(\[d\_depot,d\_central,d\_immune,d\_effector\],dtype=float\)

OBSERVEDDATAVISUALISATIONS

Thefollowingplotssummarizetheobservedtrajectories\.Examine

themcarefullybeforereadingtheODEstructureorproposingany

priors\.Usetrajectoryshape,inter\-trajectoryvariability,and

cross\-variablecorrelationstoanchoryourreasoning\.

\-mean\_trajectory\.png:mean\+/\-CIofdepot,central,immune,

andeffectorovertin\[0,50\]

\-diff\_corr\_heatmap\.png:first\-differenceSpearmancorrelations

betweenvariabletrajectorypairs

\-faceted\_by\_baseline\.png:trajectoriesstratifiedbyeachvariable’s

initialconditionquartile\(Q1lowest,

Q4highest\)

PLANNINGRULES\(CRITICAL\)

\-Examinetheobservedvisualisationplotsbeforeparsinganyparameter

\-Infereachparameter’srolestrictlyfromtheODEstructureand

theobservedtrajectorypatterns

PRIORPOLICY\(HARD\)

\-AllparametersareLogNormal;reportlinear\-spacemeanandSD\-\-\-

donotuselog\-spacemeansorSDs

\-Everyparameterentrymustinclude:mean,sd,units,andrationale

\-rationalemustbeonesentencecitingboththeinferredroleof

theparameterintheODEandamagnitudeortimescalecheck

OUTPUTCONTRACT

ReturnasingleJSONblock:

\{

"summary":"Onesentencedescribingtheoverallcalibration

strategyandkeyassumptions\.",

"params":\{

"0":\{"mean":\.\.\.,"sd":\.\.\.,"units":"\.\.\.","rationale":"\.\.\."\},

\.\.\.

\}

\}

OUTPUTCHECKLIST\(FOLLOWINORDER\)

1\.Examinethethreeobservedvisualisationplots

2\.Parseeachparameter’srolefromtheODEstructure

3\.Anchoreachpriortotheobservedtrajectorypatterns

4\.Confirmallpriorsareinlinearspace

5\.Returntheoutputcontract\-\-\-noadditionalproseafter

### G\.3Diagnosis–Update Prompts and Example

Diagnosis PromptYouareoperatingundera5\-turnbudget\.Beforeeveryresponse,

check:ifthisisturn4orlater,returntheoutputcontract

JSONimmediatelyregardlessofwhethertoolcallsarecomplete\.

ApartialanswerwithJSONisalwaysbetterthannoJSON\.

Youareanexpertdynamicalsystemsanalystwithdeepexpertisein

mechanisticmodelingandBayesianparameterestimationfor

coupledpharmacologicalandimmunologicalsystems\.

Yourtaskistoidentifythetopfailuremodesbetweentheobserved

andsynthetictrajectoriesofdepot,central,immune,andeffector

overtin\[0,50\],usingtoolcallstoconfirmeachfindingwith

computedevidence\.

Simulation:tindays,horizon=50

Statevariables:depot\(drugdepotamount\),central\(drugcentral

compartmentamount\),immune\(immunecellpopulation\),effector

\(effectorcellpopulation\)

defsystem\(biomarkers:np\.ndarray,params:np\.ndarray\)\-\>np\.ndarray:

depot,central,immune,effector=biomarkers

\#params\[0\]=absorptionratefromdepottocentral

\#params\[1\]=centraleliminationrate

\#params\[2\]=drugeffectscaleonimmunesuppression

\#params\[3\]=drughalf\-saturationconstant

\#params\[4\]=immunebaseline/sourcerate

\#params\[5\]=immunenaturaldecayrate

\#params\[6\]=effectorrecruitmentratefromimmunecontact

\#params\[7\]=effectornaturaldecayrate

k\_a=params\[0\]

k\_el=params\[1\]

e\_max=params\[2\]

ec50=params\[3\]

s\_imm=params\[4\]

d\_imm=params\[5\]

k\_recruit=params\[6\]

d\_eff=params\[7\]

drug\_effect=e\_max\*central/\(ec50\+central\+1e\-12\)

d\_depot=\-k\_a\*depot

d\_central=k\_a\*depot\-k\_el\*central

d\_immune=s\_imm\-d\_imm\*immune\-drug\_effect\*immune

d\_effector=k\_recruit\*immune\*effector\-d\_eff\*effector

returnnp\.array\(\[d\_depot,d\_central,d\_immune,d\_effector\],dtype=float\)

STATINTERPRETATIONNOTE

Allsingle\-variablestatisticsarecomputedpertrajectorythen

averagedacrosstrajectories\.Cross\-variablestatisticsare

computedpertrajectorypairthenaveraged\.

COMPARISONVISUALISATIONS\(observedvssynthetic\)

Thefollowingplotscomparetheobservedtrajectoriesagainstthe

synthetictrajectoriesgeneratedfromthebestperformingparameter

distributionssofar\(highestlogSL\)\.

\-mean\_trajectory\.png:mean\+/\-CIovertimepervariable;

bluesolid=observed,reddashed=synthetic

\-diff\_corr\_heatmap\.png:first\-differenceSpearmancorrelations;

leftpanel=observed,rightpanel=synthetic

\-faceted\_by\_baseline\.png:mean\+/\-CItrajectoriesstratifiedby

eachvariable’sinitialconditionquartile

\(Q1lowest,Q4highest\)

PARAMETERDISTRIBUTIONS

AllparametersfollowLogNormaldistributions\.Meansandstandard

deviationsarereportedinlinearspace\.

INSTRUCTIONS\(CRITICAL\)

ReasonfromtheODEstructureasyourprimarydiagnosticframework\.

Beforeexaminingtheplots,considerwhatthesystemimpliesforeach

variable:expectedtrajectoryshape,sourcesofinter\-trajectory

heterogeneity,andcouplingbetweendepot,central,immune,andeffector\.

Usethecomparisonplotstogroundthisreasoningintheobserved

trajectories\.Thencompareagainstthesynthetictrajectories:

identifywherethesyntheticdistributionfailstoreproducethe

observeddynamics,andreasonfromtheODEstructureaboutwhich

parameterrolesaremostlikelyresponsibleforeachmisalignment\.

OUTPUTCONTRACT

Usetoolcallstocomputestatisticsthatquantifysuspectedfailure

modes\.Afterreceivingresults,confirmorreviseyourassessment,

callfurthertoolsifneeded,thenreturn:

\{

"failure\_modes":\[

\{

"description":"\.\.\.",

"severity":"high\|medium\|low",

"affected\_variable":"\.\.\.",

"supporting\_evidence":"Onesentencecitingspecifictoolresults\.",

"key\_stats":\[

\{

"tool":"\.\.\.",

"variable":"\.\.\.",

"observed":\.\.\.,

"synthetic":\.\.\.,

"interpretation":"Onesentenceexplainingwhatthisgapmeans

mechanistically\."

\}

\]

\}

\]

\}

Round 1: Tool Calls and Statistics ReturnedTOOLCALLS:

mean\(depot\),mean\(central\),mean\(immune\),mean\(effector\)

std\(depot\),std\(central\),std\(immune\),std\(effector\)

level\_corr\(immune,effector\),diff\_corr\(central,immune\)

diff\_corr\(depot,central\),spearman\_trend\(immune\)

spearman\_trend\(effector\),value\_range\(immune\),value\_range\(effector\)

STATISTICSRETURNED:

statvariableobservedsyntheticgaprel\_gap

\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-

meandepot0\.1330\.211\+0\.078\+58\.4%

meancentral0\.3490\.425\+0\.076\+21\.7%

meanimmune5\.6492\.692\-2\.957\-52\.3%

meaneffector1\.6650\.681\-0\.984\-59\.1%

stddepot0\.6220\.742\+0\.120\+19\.2%

stdcentral0\.7320\.743\+0\.011\+1\.5%

stdimmune1\.2981\.800\+0\.502\+38\.7%

stdeffector0\.4570\.147\-0\.310\-67\.8%

level\_corrimmunevseffector0\.506\-0\.294\-0\.800\-158\.2%

diff\_corrcentralvsimmune\-0\.5460\.351\+0\.897\+164\.3%

diff\_corrdepotvscentral0\.5780\.705\+0\.127\+22\.0%

spearman\_trendimmune0\.8820\.800\-0\.082\-9\.2%

spearman\_trendeffector0\.518\-0\.527\-1\.045\-201\.8%

value\_rangeimmune4\.5396\.272\+1\.733\+38\.2%

value\_rangeeffector1\.4750\.482\-0\.994\-67\.4%

Round 2: Tool Calls and Statistics ReturnedTOOLCALLS:

acf\_lag1\(immune\),acf\_lag1\(effector\)

mean\_crossing\_rate\(immune\),mean\_crossing\_rate\(effector\)

turning\_point\_rate\(immune\),turning\_point\_rate\(effector\)

spectral\_entropy\(immune\),spectral\_entropy\(effector\)

STATISTICSRETURNED:

statvariableobservedsyntheticgaprel\_gap

\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-

acf\_lag1immune0\.9630\.950\-0\.014\-1\.4%

acf\_lag1effector0\.9670\.972\+0\.005\+0\.5%

mean\_crossing\_rateimmune0\.0180\.020\+0\.003\+14\.3%

mean\_crossing\_rateeffector0\.0130\.020\+0\.008\+60\.0%

turning\_point\_rateimmune0\.0110\.010\-0\.001\-5\.7%

turning\_point\_rateeffector0\.0100\.020\+0\.010\+104\.1%

spectral\_entropyimmune1\.2090\.912\-0\.297\-24\.6%

spectral\_entropyeffector1\.2610\.665\-0\.596\-47\.3%

Round 3: Tool Calls and Statistics ReturnedTOOLCALLS:

dfa\_alpha\(immune\),dfa\_alpha\(effector\)

spearman\_trend\(depot\),spearman\_trend\(central\)

level\_corr\(central,immune\),level\_corr\(central,effector\)

diff\_corr\(immune,effector\),diff\_corr\(depot,central\)

STATISTICSRETURNED:

statvariableobservedsyntheticgaprel\_gap

\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-

dfa\_alphaimmune1\.5661\.408\-0\.158\-10\.1%

dfa\_alphaeffector2\.0592\.086\+0\.026\+1\.3%

spearman\_trenddepot\-0\.841\-1\.000\-0\.159\-19\.0%

spearman\_trendcentral\-0\.941\-0\.939\+0\.002\+0\.2%

level\_corrcentralvsimmune\-0\.912\-0\.878\+0\.034\+3\.7%

level\_corrcentralvseffector\-0\.5030\.472\+0\.975\+193\.7%

diff\_corrimmunevseffector0\.1360\.502\+0\.366\+269\.5%

diff\_corrdepotvscentral0\.5780\.705\+0\.127\+22\.0%

Diagnosis Output\{

"failure\_modes":\[

\{

"description":"Effectordynamicsarequalitativelywrong:synthetic

trajectoriescollapsetoamuchlower,lessvariable,andoppositely

trendingeffectorpopulation,indicatingthemodelfailstocapture

theimmune\-to\-effectorrecruitmentfeedback\.",

"severity":"high",

"affected\_variable":"effector",

"supporting\_evidence":"Effectormean1\.665vs0\.681,std0\.457vs

0\.147,Spearmantrendflipsfrom\+0\.518to\-0\.527\.",

"key\_stats":\[

\{"tool":"mean","variable":"effector",

"observed":1\.665,"synthetic":0\.681,

"interpretation":"Syntheticeffectorburdenfartoolow,consistent

withinsufficientrecruitmentoroverlystrongdecay\."\},

\{"tool":"spearman\_trend","variable":"effector",

"observed":0\.518,"synthetic":\-0\.527,

"interpretation":"Signreversalimpliessyntheticsystemdrives

effectordeclinewhenobservedsystemgenerallyrises\."\}

\]

\},

\{

"description":"Immune\-effectorcouplingmis\-specified:positive

co\-movementobservedbutanti\-correlationinsynthetictrajectories\.",

"severity":"high",

"affected\_variable":"immune",

"supporting\_evidence":"Levelcorrelation\+0\.506observedvs\-0\.294

synthetic;diffcorrelation0\.136vs0\.502\.",

"key\_stats":\[

\{"tool":"level\_corr","variable":"immunevseffector",

"observed":0\.506,"synthetic":\-0\.294,

"interpretation":"Syntheticmodelreversessteady\-statealignment

betweenimmuneandeffector\."\},

\{"tool":"spectral\_entropy","variable":"immune",

"observed":1\.209,"synthetic":0\.912,

"interpretation":"Syntheticimmuneseriestoospectrallysimple,

consistentwithmissingheterogeneity\."\}

\]

\},

\{

"description":"Immunedynamicstoosuppressedandover\-dispersed:

meanfarbelowobserved,variabilityinflated\.",

"severity":"high",

"affected\_variable":"immune",

"supporting\_evidence":"Immunemean5\.649vs2\.692;std1\.298vs1\.800;

spectralentropy1\.209vs0\.912\.",

"key\_stats":\[

\{"tool":"mean","variable":"immune",

"observed":5\.649,"synthetic":2\.692,

"interpretation":"Syntheticimmuneabundancetoolow,pointingto

excessivesuppressionorinsufficientreplenishment\."\}

\]

\},

\{

"description":"Depot\-centralPKphasedirectionallycorrectbut

quantitativelytoohighinabsolutemass\.",

"severity":"medium",

"affected\_variable":"depot",

"key\_stats":\[

\{"tool":"mean","variable":"depot",

"observed":0\.133,"synthetic":0\.211,

"interpretation":"Syntheticdepotpersistencetoostrong,

consistentwithabsorptionbeingtooslow\."\}

\]

\},

\{

"description":"Centralbroadlycapturedbutdownstreamcouplingto

effectorhaswrongsign\.",

"severity":"medium",

"affected\_variable":"central",

"key\_stats":\[

\{"tool":"level\_corr","variable\_x":"central",

"variable\_y":"effector",

"observed":\-0\.503,"synthetic":0\.472,

"interpretation":"Modelassignswrongsignofexposure\-response

association,amechanisticredflagforthePDcoupling\."\}

\]

\}

\]

\}

Update PromptYouareoperatingundera5\-turnbudget\.Beforeeveryresponse,

check:ifthisisturn4orlater,returntheoutputcontract

JSONimmediatelyregardlessofwhethertoolcallsarecomplete\.

ApartialanswerwithJSONisalwaysbetterthannoJSON\.

Youareanexpertdynamicalsystemsanalystwithdeepexpertisein

mechanisticmodelingandparameterdistributioninferencefor

coupledpharmacologicalandimmunologicalsystems\.

YourtaskistoproposeupdatedLogNormalparameterdistributions

thatrecoverthetrueunderlyingdynamicsoftheobserved

trajectories\.ParameterchangesmustbeconsistentwiththeODE

structureandsupportedbyevidencefromiterationhistory\.

Simulation:tindays,horizon=50

Statevariables:depot\(drugdepotamount\),central\(drugcentral

compartmentamount\),immune\(immunecellpopulation\),effector

\(effectorcellpopulation\)

defsystem\(biomarkers:np\.ndarray,params:np\.ndarray\)\-\>np\.ndarray:

depot,central,immune,effector=biomarkers

\#params\[0\]=absorptionratefromdepottocentral

\#params\[1\]=centraleliminationrate

\#params\[2\]=drugeffectscaleonimmunesuppression

\#params\[3\]=drughalf\-saturationconstant

\#params\[4\]=immunebaseline/sourcerate

\#params\[5\]=immunenaturaldecayrate

\#params\[6\]=effectorrecruitmentratefromimmunecontact

\#params\[7\]=effectornaturaldecayrate

k\_a=params\[0\]

k\_el=params\[1\]

e\_max=params\[2\]

ec50=params\[3\]

s\_imm=params\[4\]

d\_imm=params\[5\]

k\_recruit=params\[6\]

d\_eff=params\[7\]

drug\_effect=e\_max\*central/\(ec50\+central\+1e\-12\)

d\_depot=\-k\_a\*depot

d\_central=k\_a\*depot\-k\_el\*central

d\_immune=s\_imm\-d\_imm\*immune\-drug\_effect\*immune

d\_effector=k\_recruit\*immune\*effector\-d\_eff\*effector

returnnp\.array\(\[d\_depot,d\_central,d\_immune,d\_effector\],dtype=float\)

CALIBRATIONGOAL\(CRITICAL\)

Addressdiagnosedfailuremodesstrictlyinorderofseverity\-\-\-

highbeforemediumbeforelow\.Donotmakechangesthatimprovea

lower\-severitymismatchatthecostofworseningahigher\-severity

one\.Reasonatthelevelofparameterrolesbeforereasoningabout

specificvalues\-\-\-whateachparameterdoesinthesystemmatters

morethanitsindex\.

DIAGNOSEDFAILUREMODES

Thefollowingfailuremodesandsummarystatisticswerediagnosed

fromthebestperformingparameterdistributionswhichhavethe

highestlogSLvaluesofar\.

\[\.\.\.failuremodesJSONfromdiagnosisoutputabove\.\.\.\]

MANDATORYWORKFLOW

STEP1\-\-\-Mapfailuremodestoparameterroles

Foreachdiagnosedfailuremode:

\-identifywhichdynamicalpropertyoftheODEismisrepresented

\-identifywhichparameterrolesintheODEgovernthatproperty

\-explainhowthoserolesinteractjointly\-\-\-adjustingonerole

willaffectotherpropertiesthroughsharedODEterms

\-statethedirectionofcoordinatedadjustmentneededacross

allgoverningrolessimultaneously

Reasonatthelevelofparameterrolesbeforereasoningabout

specificvalues\.

STEP2\-\-\-Checkforroleconflicts

Identifywhethercorrectingonefailuremodethroughitsgoverning

parameterroleswillworsenanotherfailuremodethroughsharedor

coupledroles\.Ifaconflictexists,prioritizethehighest

severityfailuremodeandacceptthetradeoffexplicitly\.

STEP3\-\-\-Reviewiterationhistory

Thebest\-performingiterationsareprovidedinthecontextmessage

above\.ForeachroleadjustmentidentifiedinSteps1and2,find

whathappenedwhensimilarroleadjustmentsweremadepreviously\.

Foreveryroleadjustmentyouplantomake,cite:oneiteration

whereasimilaradjustmenthelped,onewhereithurtorcauseda

conflict,andtheODE\-basedreasonwhy\.Ifthehistoryistoothin

toestablishthis,preferthesmallestpossiblechange\.

STEP4\-\-\-Proposecoherentjointupdate

Proposeparametervaluesthatjointlyaddressthehighestpriority

failuremodes\.Everychangemust:

\-correctaspecificdynamicalpropertythroughitsgoverningroles

\-beconsistentwithhowthoserolesinteractintheODE

\-becoordinatedacrossallparametersservingrelatedroles\-\-\-

notone\-at\-a\-timecorrections

\-beminimal\-\-\-preferthesmallestchangethataddressesthe

failuremode,sincelargechangesinnonlinearsystemshave

unpredictablejointeffects

STABILITYCHECK\(HARD\)

Beforefinalising,verifythattheproposedparametervalues

supportastablesystemequilibrium\.Iftheproposedparameters

woulddrivethesystemtocollapseorexplosion,revisethem

beforereturningtheoutputcontract\.

ITERATIONHISTORY

Anindexofallpastiterationsisshownbelow\.Thebest\-performing

iterations\(fullparamsandlogSLbreakdown\)areprovidedinthe

contextmessageabove\-\-\-usethemtogroundyourparameterchanges\.

itertypelog\_slsummary

0initial\_inferenceN/AModeratelyinformativepriors

centeredonfastdepotabsorption,

moderatecentralelimination,

saturablebutnotoverwhelmingdrug

suppressionofimmunecells,and

slow\-to\-moderateeffectorexpansion,

chosentoreproducetheobserved

rapiddrugwashout,U\-shapedimmune

rebound,andgradualeffectorrise

over50days\.

1constraint\_violationN/ACorrectedthedepotabsorption/

eliminationandeffectorgrowth/

decaypriors,becausetheinvalid

trajectoriesaredominatedbytiny

negativedepotvaluesfromoverly

aggressivedepletionandextreme

effectorexplosionsfromtoo\-strong

netrecruitment\.

2normal\-1865719\.7Correctedpriorsfordepot

absorptionandeffectorgrowth/decay

becausetheinvalidtrajectoriesare

drivenbyslightdepotnegativity

andrunawayeffectoramplification,

whiletheremainingparameterswere

leftunchangedastheyarebroadly

consistentwiththeobservedranges\.

3normal\-7530\.6Reducednetimmunosuppressionand

strengthenedeffectorpersistence

whileslightlyeasingdepotexposure

,becausethedominanterrorisan

immune/effectorcollapsewithsign\-

reversedcouplingthatmustbefixed

beforeanysecondaryPKsmoothing

adjustments\.

PRIORPOLICY\(HARD\)

\-AllparametersareLogNormal;reportlinear\-spacemeanandSD\-\-\-

donotuselog\-spacemeansorSDs

\-Everyparameterentrymustinclude:mean,sd,units,andrationale

\-rationalemustbeonesentencecitingtheparameterroleinthe

ODEandthecoordinatedcorrectionmade\-\-\-nottheparameterindex

OUTPUTCONTRACT

ReturnasingleJSONblockcontainingALLparameters:

\{

"summary":"Onesentencedescribingthekeyroleadjustments

madeandwhy\.",

"params":\{

"0":\{"mean":\.\.\.,"sd":\.\.\.,"units":"\.\.\.","rationale":"\.\.\."\},

\.\.\.

\}

\}

Update Output\{

"summary":"Increaseimmunereplenishmentandrelaxdrug\-mediated

suppressionwhilestrengtheningeffectorpersistenceviaconservative

recruitment/decaybalance,toaddresstheunderpoweredsign\-flipped

effectorresponseandsuppressedimmunelevel\.",

"params":\{

"0":\{"mean":0\.500,"sd":0\.100,"units":"1/day",

"rationale":"Slightlyreducesdepotpersistencetoaddressmodest

PKoverpredictionwithoutdisturbingthePKskeleton\."\},

"1":\{"mean":0\.240,"sd":0\.050,"units":"1/day",

"rationale":"Unchanged;centraleliminationbroadlyconsistent

withobservedtrendandvariability\."\},

"2":\{"mean":0\.900,"sd":0\.300,"units":"dimensionless",

"rationale":"Reducedtoweakenimmunesuppressionind\_immune,

raisingimmunelevelsandrestoringcorrectcoupling\."\},

"3":\{"mean":0\.800,"sd":0\.200,"units":"centralunits",

"rationale":"Increasedtodelaysaturationofdrugeffect,

coordinatingwithlowereffectscale\."\},

"4":\{"mean":0\.340,"sd":0\.100,"units":"immuneunits/day",

"rationale":"Raisedtoliftimmunebaselinetowardobservedmean,

counteractingexcessivedepletion\."\},

"5":\{"mean":0\.030,"sd":0\.010,"units":"1/day",

"rationale":"Slightlyreducedtosupporthighersustainedimmune

levelwithoutdestabilisingthesystem\."\},

"6":\{"mean":0\.012,"sd":0\.004,"units":"1/\(immune\*day\)",

"rationale":"Increasedtostrengthenimmune\-to\-effectorcoupling

soeffectorriseswithimmuneavailability\."\},

"7":\{"mean":0\.032,"sd":0\.010,"units":"1/day",

"rationale":"Loweredtoimproveeffectorpersistenceandcorrect

thesign\-reverseddownwardtrend\."\}

\}

\}

### G\.4Contamination Analysis Prompt

Contamination Analysis Prompt \(PKPD\-Immune\)Youareanexpertindynamicalsystemsmodeling\.

Observablestatevariables:

\-depot\(drugamount,a\.u\.\):druginabsorptioncompartment

\-central\(drugamount,a\.u\.\):drugincentralcompartment

\-immune\(cellcount,a\.u\.\):immunecellpopulation

\-effector\(cellcount,a\.u\.\):effectorcellpopulation

Timeunit:days\.Observationsspantin\[0,50\]\.

Task:ProposethemostlikelyODEstructuregoverning

d\_depot/dt,d\_central/dt,d\_immune/dtand

d\_effector/dt\.Writetheequationsexplicitly

andidentifyallparameters\.