Data-Driven Discovery of Dynamical Systems in Pharmacology using Large Language Models

Authors: samuel  holt
zhaozhi  qian
james  weatherall
astrazeneca  mihaela van der schaar
resp  llm
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-14 14:42:51

Abstract: The discovery of dynamical systems is crucial across a range of fields, including pharmacology, epidemiology, and physical sciences. *Accurate* and *interpretable* modeling of these systems is essential for understanding complex temporal processes, optimizing interventions, and minimizing adverse effects. In pharmacology, for example, precise modeling of drug dynamics is vital to maximize therapeutic efficacy while minimizing patient harm, as in chemotherapy. However, current models, often developed by human experts, are limited by high cost, lack of scalability, and restriction to existing human knowledge. In this paper, we present the **Data-Driven Discovery (D3)** framework, a novel approach leveraging Large Language Models (LLMs) to iteratively discover and refine interpretable models of dynamical systems, demonstrated here with pharmacological applications. Unlike traditional methods, D3 enables the LLM to propose, acquire, and integrate new features, validate, and compare dynamical systems models, uncovering new insights into pharmacokinetics. Experiments on a pharmacokinetic Warfarin dataset reveal that D3 identifies a new plausible model that is well-fitting, highlighting its potential for precision dosing in clinical applications.

Main paper: 
 introduction: 
 The discovery of dynamical systems models plays a fundamental role across various domains, including pharmacology, epidemiology, and physical systems. In pharmacology, pharmacokinetic (PK) models are essential for understanding and predicting the time course of drug absorption, distribution, metabolism, and excretion in the body, which is crucial for optimizing therapeutic efficacy, minimizing toxicity, and personalized treatment regimens in diseases such as cancer, cardiovascular disorders, and infectious diseases [Gabrielsson and Weiner, 2001, Hedaya, 2012] . For example, cardiovascular disorders requiring Warfarin therapy affect tens of millions of individuals globally [Lee and Klein, 2013] , highlighting the critical need for superior models to understand the dynamic impact of drugs and precision dosing, thereby assisting both doctors and patients.The goal of PK modeling is to discover an underlying closed-form ordinary differential equation (ODE) f from a dataset of observed patient trajectories. This problem is of significant interest to the machine learning (ML) community, as evidenced by previous non-interpretable ML modeling approaches aimed at developing better PK models [Chan and van der Schaar, 2022] . Such datasets are collected from expensive human clinical trials, necessitating the careful selection of which features to acquire and analyze during these trials [Guidance, 2010] . Pharmacometricians then leverage this data, using their existing knowledge to select appropriate pharmacokinetic models and employ standard statistical tools to infer the parameters of these models. For instance, they often fit a single compartmental PK model, a closed-form ODE model with 3-5 parameters, to the data [Chen and Abuassba, 2021] . This process is iterative, refining the models until the most accurate one that generalizes well to new patient trajectories is identified. However, this traditional model discovery approach is often ad-hoc and fundamentally limited by the human expert's time, experience, and implicit assumptions about the underlying pharmacokinetic processes.Discovering interpretable pharmacokinetic models is traditionally performed by human experts through a scientific process of proposing, collecting, and validating models. Recently, there has been a growing call for artificial intelligence-driven methods to enhance pharmacological modeling [Ryan et al., 2024 , Singh et al., 2023 , Cheng et al., 2022] . This established process involves iterative steps, each presenting distinct challenges in constructing accurate pharmacokinetic dynamics models.Selecting appropriate temporal models involves complex decisions on model class, state and feature variables, parameters, and their relationships, far exceeding the simplicity of standard pharmacometrics models like single, double, or triple compartmental models [Chen and Abuassba, 2021] . This complexity demands aligning model intricacy with data availability, where simple models suit small datasets and complex models fit larger ones. Acquiring the right features is crucial, as it must be done with limited prior information, unlike feature selection, which evaluates pre-existing features [Li et al., 2017] . The evaluation must identify inadequate models and explain why, requiring diverse and timeintensive tools beyond a single metric like validation MSE. These challenges are interrelated: strategic data acquisition and thorough evaluation are essential for effective modeling. Despite advancements in tackling these challenges independently, temporal modeling remains manual, limiting speed and scalability, underscoring the need for integrated and automated systems to enhance efficiency.To address these challenges, we develop the Data-Driven Discovery (D3) framework, leveraging Large Language Models (LLMs) to iteratively discover and refine interpretable dynamics models, directly relevant to pharmacology, epidemiology and ecology applications. D3 features three agents: Modeling, Feature Acquisition, and Evaluation, which collaborate iteratively. The Modeling Agent uses LLMs' code generation and natural language understanding to explore the model space by generating hypothetical models. The Data Acquisition Agent utilizes LLMs' zero-shot and few-shot learning to optimize data acquisition based on summary statistics, text descriptions, and prior knowledge. The Evaluation Agent conducts comprehensive evaluations by computing refined validation MSE metrics, and alignment with prior process understanding, providing precise feedback to improve the other agents' performance.⃝ We propose the Data-Driven Discovery (D3) framework, a novel approach leveraging Large Language Models (LLMs) to iteratively discover and refine interpretable dynamics models, advancing pharmacokinetic modeling (Section 3). 2 ⃝ D3 overcomes the challenges of uncovering interpretable dynamical systems by using LLMs to explore vast model spaces and integrate unstructured data, producing models with few parameters that rival black box neural networks. D3 accurately discovers dynamics models across pharmacology, epidemiology, and ecology, matching the accuracy of existing methods while enhancing interpretability. It discovers a new, more accurate PK model for Warfarin, validated by expert pharmacometricians. We also gain insight into D3's ability to iteratively improve its models, selectively acquire features, and discover precise dynamics models.Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?Answer: [Yes] Justification: Yes, in this paper we introduce the data-driven discovery (D3) framework, verify it experimentally against existing modeling ML methods, and verify the new novel discovered PK model on Wafarin with pharmacometrician experts. We detail the framework in Section 3, provide full implementation details in Appendix F, provide the experimental results to verify our core claims of discovering well-fitting models in Section 6, and our discovered Wafarin PK model as a case study in Section 6.1. We further discuss key assumptions of the method in Section 3 and outline important limitations in Section 7, with the "limitations & future work section".Guidelines:• The answer NA means that the abstract and introduction do not include the claims made in the paper. • The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or NA answer to this question will not be perceived well by the reviewers. • The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings. • It is fine to include aspirational goals as motivation as long as it is clear that these goals are not attained by the paper. 
the vast model space for temporal modeling: 
 In this section, we introduce the model space that D3 searches through, which contains a diverse range of temporal models with varying degrees of refinement. Our focus is on Ordinary Differential Equation (ODE) models, which is one of the most widely used methods for modeling temporal dynamics, which include pharmacokinetic, physiological and epidemiological dynamics [Auger et al., 2008] . We are interested in modeling the evolution of a set of state variables x n (t) ∈ R D for individuals n = 1, . . . , N over the time horizon t ∈ [0, T ]. The state variables represent the target variables that the human experts are interested in modeling; and there may exist other feature variables that can help predict the evolution of state variables. We assume we are provided with a dataset of individuals trajectories such that D = {(x n (t), a n (t)) | n = 1, . . . , N, t ∈ [0, T ]}, where a n (t) ∈ R K denotes the observed individual-level features for individual. We can categorize ODE models into three levels (R1 -R3) based on their capacity for refinement.R1: Refinement through initial conditions x n (0). This level involves the simplest form of ODEs, commonly referenced in scientific literature [Schiesser, 2014] . These equations are defined over state variables x(t) ∈ R D , incorporate global parameters θ, and involve a function f : R D+1 → R D . The time t spans the interval [0, T ]: ẋn (t) = f (x n (t), t, θ), ∀n ∈ [N ] (1)Considering N individuals, each represented by the trajectory x n (t) for n = 1, . . . , N , all follow the equation above. According to the uniqueness theorem for initial value problems in ODEs [Lindelöf, 1894] , any differences among individuals' trajectories can be attributed solely to variations in the initial conditions x n (0), under common regularity conditions on f . R2: Refinement through observed features a n . Suppose we have access to a set of individual-level feature variables a n . We can enhance the model's refinement by incorporating these features into the ODE to enable unique dynamics for different individuals:EQUATIONR3: Refinement through acquired features h n . Often, not all relevant variables are initially measured and available for analysis. In such cases, one may consider acquiring additional features h n (t) ∈ R J to be integrated into the ODE. Here h n (t) is a collection of J features from the set of all acquirable features H, where J ≤ |H|. This third level of refinement reflects one data acquisition challenge, i.e., identifying and collecting informative features h n (t):EQUATIONAfter selecting the level of refinement (R1-R3), it is necessary to determine the functional form of the ODE f , which defines the interactions among various variables and parameters that govern the dynamics. There are two dominant approaches to parameterize f : using a concise closed-form white-box equation f ∈ M C or employing a neural network f ∈ M N . The former approach is prevalent in the field of symbolic regression [Billard and Diday, 2002] , while the latter is utilized in Neural ODEs [Chen et al., 2018] . Additionally, f can incorporate both a closed-form component and a neural component, resulting in a hybrid ODE model. Nevertheless, all three approaches involve searching within large combinatorial spaces-closed-form equations, neural architectures, and both-presenting significant computational challenges.3 Data-Driven Discovery (D3)We aim to uncover the true underlying pharmacokinetic dynamical system, denoted as f * . Our goal is to continually approximate f * as accurately as possible using proposed models f θ , with parameters θ. For simplicity, we will refer to these models as f in all subsequent references, omitting the θ term.To appropriately identify the refinement model level needed from the vast model space discussed above, we now introduce the Data-Driven Discovery (D3) framework, as depicted in Figure 1 . This consists of three LLM agents: the Modeling Agent G, the Feature Acquisition Agent A, and the Evaluation Agent E. These three agents work together to form the model improvement loop and the data acquisition loop. The implementation details of the Agents are provided in Appendix F.Question: Does the paper discuss the limitations of the work performed by the authors?Answer: [Yes] Justification: We have a dedicated limitations section in Section 7, titled "limitations & future work section". We also discuss the key assumptions of the method when we introduce it, such as in Section 3, for example in Section 3.3 we state the assumptions made, and how relevant they are in clinical practice, where D3 is designed to work. We also discuss the computational efficiency of D3 in Appendix F.Guidelines:• The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper. • The authors are encouraged to create a separate "Limitations" section in their paper.• The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be. • The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on a few datasets or with a few runs. In general, empirical results often depend on implicit assumptions, which should be articulated. • The authors should reflect on the factors that influence the performance of the approach.For example, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon. • The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size. • If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness.• While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations. 
inputs to data-driven discovery (d3): 
 To begin the process, D3 requires a clearly defined description of the system in natural language, that details the modeling task, of creating a well-fitting model of either a white-box model or a hybrid model; that is a white-box model with a black-box neural network fitted to the residuals of the white-box model. D3 supports both discovering only white-box models and hybrid models, and we leave this choice up to the end user to decide for their application. Specifically, a user must provide a dataset of individual's trajectories D, and a system description and the names and descriptions of any features that exist within the dataset, which can include the ranges of those features. Prior information can be added here by the practitioner into the system description if it is available. We now discuss how the three agents of Modeling, Feature Acquisition and Evaluation interact to discover iteratively better-fitting f models. We provide full implementation details for the framework in Appendix F.D3's Modeling Agent G is tasked to iteratively propose and refine the temporal model f i , where i = 1, . . . , I indicates the iteration. We design G following the three principles below. D3 can discover a well-fitting model f , that is either a white-box model or a hybrid model, combining a white-box model with a neural network component fit to the residuals. First, a user provides a system description c for the model, which includes the feature names, their units and ranges. The Modeling agent uses the textual description and generates a model, represented as code. Next, the Evaluator agent evaluates the generated model on a held-out validation dataset and provides a loss metric. Next, the Feature Acquisition agent decides if it wants to acquire a new feature, that could aid in generating better models. This iterative process repeats, when there exists one or more previously generated models, the Evaluator agent provides a verbal reflection on how the model can be improved, which is used by the Modeling agent at subsequent iterations. This process repeats, discovering better fitting models, and after a set number of iterations, the best fitting model is returned, as code with its parameters optimized to the given training subset of the given dataset D.Representing the model f i as code. As we discussed in Section 1, one key challenge in evidencebased refinement is the vast space of possible ODE models, ranging from closed-form ODEs to complex neural Neural ODEs and hybrid models. As a prerequisite to finding the most suitable model, we need a consistent and flexible way to represent these diverse models such that they can be trained and validated on data. As such, D3 represents the model f i as code in Python language. Specifically, each model is defined as a class derived from 'nn.Module' in PyTorch [Paszke et al., 2017] that specifies the free parameters and the computation logic. The Modeling Agent G leverages an LLM to generate the model code f i in each step.Informing generation with context c. Effectively searching through the vast model space presents another challenge. Thankfully, in many applications there exists rich context c about prior works and domain knowledge. However, this rich source of information was previously locked in unstructured documents and texts, making it hard for automated algorithms to make use of it. By leveraging the LLM's capability to understand unstructured documents, the Modeling Agent G can now access and leverage this information, therefore proposing more informed models and searching more effectively. In this work, we consider the provision of the following categories of problem context of {system description} and {feature description} (Appendix F).Iterative model improvements based on rich verbal feedback r i-1 and memory s i-1 . To enhance the quality of the generated model f i , we have designed the agent G to iteratively improve the model based on the presence of any evaluation verbal feedback and its knowledge of the previously generated models and how they performed on a validation subset of given dataset D of trajectories. At each iteration i, the Evaluation Agent E evaluates the previous model f i-1 using the validation data subset of D i and provides rich verbal feedback r i-1 (c.f. Figure 1 ), which forms the foundation for the next generated model f i . The Modelling Agent G utilizes memory s i-1 to track the top-k best-performing models so far, along with their associated evaluation feedback, here their respective losses. Further details on evaluation and data acquisition are discussed in the following section. This iterative process allows for continuous refinement and optimization of models based on their past performance and targeted feedback.In summary, the Modeling Agent G encapsulates a comprehensive range of model structures through computer code f i and efficiently explores the vast model space by integrating unstructured prior knowledge c, and utilizing data-informed feedback r i-1 along with memory s i-1 . The operation of G can be mathematically represented as follows:f i = G(c, f i-1 , D i , r i-1 , s i-1 )(4)D3's Feature Acquisition Agent A is tasked with iteratively proposing new features for acquisition. These features are selected based on their potential to improve the next model's performance in modeling the target state x n (t).EQUATIONWhere h i is the next feature to acquire. Specifically, when a feature is acquired we make the assumption, that it is acquired for all individuals, therefore at the next iteration the dimension of the individual-level features a n (t) ∈ R K , n = 1, . . . , N increases by one, to a n (t) ∈ R K+1 , n = 1, . . . , N . It is possible to do this in clinical trials, where pharmacokinetic datasets of patients are collected [Dziura et al., 2013] , for example acquiring the individual's ages, or other lab biomarker measurements or tests [Derraik et al., 2021] .Estimating the value of the feature. We leverage the value of information [Feltham, 1968] framework to provide a principled way to determine which feature, if any to acquire next. This follows as:EQUATIONwhere the first two terms capture the improvements in validation loss L when model f i is trained on the new dataset that includes the new feature, that is D = {(x n (t), a n (t), h n (t))|∀n ∈ [N ]} and l(h i ) represents the cost for acquiring h i . Estimating V (h i ) presents a statistical challenge and a computational challenge. Firstly, in practice, we often do not have access to h i to begin with (as the goal is to acquire new features)-hence we need a way to estimate E[L(f i , h i , D i-1 )] without fully accessing h i . Secondly, even when h i is available, computing the two loss terms involves fitting models twice which can be computationally challenging.To address these challenges, the Feature Acquisition Agent A predicts the value of feature V (h i ) from the available information about h i , as follows:EQUATIONwhere c hi represents the unstructured prior information about the feature h i (e.g. data descriptions and metadata), T (h i ) represents the available summary statistics (e.g. range and type of values from a small subset of data h i ), D i-1 represents the existing data, and g is the prediction function.The prediction challenge above corresponds to zero-shot or few-shot learning in ML. In the zero-shot case, we need to predict V (h i ) without having access to any "label" V (h), for some feature h. In the few-shot case, we have access to a small number of feature-label pairs, {c hi , T (h i ), h i-1 , V (h i-1 )} to learn the prediction function g. The Feature Acquisition Agent A leverages LLM's capability in zero-shot and few-shot learning to address this estimation challenge. Specifically, we inform the LLM by providing high-level statistics T (h i ) and a description of the feature c hi in the form of the feature name to warm-start the zero-shot learning. When the h i has been included in the data, we inform the LLM about the validation loss associated with those data v(h i ) and add it to the few-shot examples to guide the estimation.Context, feedback and memory. Similar to G, A is also able to leverage unstructured problem context c, evaluation feedback r i-1 , and memory s i-1 to guide the data acquisition proposal.The Evaluation Agent E is vital for iteratively improving the Modeling Agent G and Feature Acquisition Agent A. It provides feedback r i :EQUATIONFeedback r i can be numeric, such as validation loss or mean squared error (MSE) of the model f i on the validation subset D i , or more refined metrics like validation loss per target dimension in D. It can also be textual, offering detailed suggestions for model improvements, including clinically plausible modifications generated by an LLM [Shinn et al., 2024] .Evaluating dynamical systems is complex due to intricate interactions and temporal dependencies. The Evaluation Agent dynamically assesses both model performance and plausibility, ensuring models are both accurate and interpretable. This comprehensive evaluation is crucial for developing sophisticated, clinically relevant pharmacokinetic models that address real-world complexities. Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof?Answer: [NA] Justification: The paper does not include any theoretical results, therefore this question is not applicable.Guidelines:• The answer NA means that the paper does not include theoretical results.• All the theorems, formulas, and proofs in the paper should be numbered and crossreferenced. • All assumptions should be clearly stated or referenced in the statement of any theorems.• The proofs can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide a short proof sketch to provide intuition. • Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material. • Theorems and Lemmas that the proof relies upon should be properly referenced. 
related works: 
 C × Neural ODE R1, R2 M N × AI Feynman R1 M C ✓ Eureka R1, R2 M C ✓ AFA × ✓ Train Active Learning × × Train Active Testing × × Eval D3 R1 -R3 M C ∪ M N ✓ ✓ ✓ BothOur work focuses on autonomously learning temporal models while acquiring data, with several relevant research strands, as summarized in Table 1 , which is expanded in Appendix A.ODE learning methods. Symbolic regression methods like SINDy and D-CODE [Koza, 1994 , Brunton et al., 2016 , Qian et al., 2022] can discover closed-form ODEs f ∈ M C using genetic algorithms. They can incorporate features a n (R2 refinement) but struggle with many variables (e.g., >20) due to computational complexity. Neural ODEs use neural networks f ∈ M N to handle many variables [Chen et al., 2018 , Dupont et al., 2019 , Zaytar and El Amrani, 2016 , Devlin et al., 2018 , Sehovac and Grolinger, 2020] , but have numerous free parameters, risking overfitting with insufficient data. D3 captures f ∈ M C ∪ M N and supports R1 -R3 refinement, adapting based on data availability. Both Symbolic regression and Neural ODEs require manual context incorporation c via hyperparameters and do not address Data Acquisition Challenges.AI for automated modeling. AI Feynman [Udrescu and Tegmark, 2020] enhances symbolic regression by using the physical units of variables, which provide additional constraints to narrow the search space. However, its utility outside physics is limited, as units in other fields carry less information.Eureqa [Ma et al., 2023] applies LLMs to model the reward function in reinforcement learning, integrating unstructured context and representing models as code. Unlike D3, Eureqa cannot automatically refine models, acquire new data, or easily apply to learning temporal dynamics like ODEs.Data Acquisition. Methods like Active Feature Acquisition (AFA) [Ma et al., 2018 , Gong et al., 2019] measure additional features on existing samples to improve performance. Active Learning [Sebastiani and Wynn, 2000 , Settles, 2009 , Sener and Savarese, 2017 , Imberg et al., 2020] and Active Testing [Lowell et al., 2018 , Kossen et al., 2021] acquire new samples to increase training or evaluation data size while keeping features constant. These methods do not address the zero-shot setting, where data acquisition decisions must be made before observing any target data (Equation 7). Additionally, most Active Learning methods focus on acquiring labels for supervised learning, which is not directly applicable to temporal modeling tasks.Question: Does the paper fully disclose all the information needed to reproduce the main experimental results of the paper to the extent that it affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?Answer: [Yes] Justification: To ensure thorough reproducibility, in addition to D3's description within the main paper (Section 3), we include full implementation details of D3 in Appendix F. To fully reproduce the experiments we detail the exact datasets, and how they were preprocessed or generated in Appendix B, how the respective evaluation metrics were computed in Appendix C, along with the corresponding full implementation details of the implemented benchmark baseline methods in Appendix E. 
experiments and evaluation: 
 In this section, we demonstrate that D3 can discover well-fitting pharmacokinetic dynamical system models for a range of diverse PK datasets, including an epidemiological and ecological dataset.Benchmark Datasets. Our evaluation encompasses six real-world datasets with clinical relevance, each originating from either real-world data or highly accurate simulators developed by human experts.Three datasets are based on a state-of-the-art biomedical Pharmacokinetic-Pharmacodynamic (PKPD) model of lung cancer tumor growth, which simulates the combined effects of chemotherapy and radiotherapy in lung cancer [Geng et al., 2017] (Equation ( 11)). This model has been widely utilized in previous research [Bica et al., 2020 , Seedat et al., 2022 , Melnychuk et al., 2022] . Specifically, we employ this bio-mathematical model to generate three variations: lung cancer without treatment (Lung Cancer), lung cancer treated with chemotherapy (Lung Cancer (with Chemo.)), and lung cancer treated with both chemotherapy and radiotherapy (Lung Cancer (with Chemo. & Radio.)).Additionally, we utilize an intricate COVID-19 epidemic agent-based simulator (COVID-19) [Kerr et al., 2021] , to provide an epidemiological dataset. Another dataset comes from an ecological model simulating a microcosm of algae, flagellate, and rotifer populations (Plankton Microcosm), replicating a three-species prey-predator experimental system [Hiltunen et al., 2013] . Finally, we include a real Pharmacokinetic (PK) dataset of Warfarin patients (Warfarin) [Janssen et al., 2022] .Detailed information about all benchmark datasets is provided in Appendix B.Benchmark Methods. To evaluate the performance of D3 and establish its competitive performance, we conduct comparisons with leading modeling methods for ODEs. Specifically, we benchmark against advanced black-box models that have many parameters, such as neural ODEs with action inputs, known as DyNODE [Chen et al., 2018 , Alvarez et al., 2020] . Also in this class are recurrent neural networks (RNN) and a state-of-the-art transformer model (Transformer). In addition, we include white-box transparent dynamical systems models identified through equation discovery techniques, such as Sparse Identification of Nonlinear Dynamics (SINDy) [Brunton et al., 2016] . Moreover, D3 supports two discovery modes, of which we compare against both, discovering only white-box models (D3-white-box) and discovering hybrid models (D3-hybrid) 3 . We also perform ablations of D3, of zero-shot generated model from D3 as (ZeroShot) and the same model with optimized parameters (ZeroOptim). Detailed descriptions of the implementations, hyperparameters, and experimental procedures for these benchmarks can be found in Appendix E.Evaluation Metrics. To assess the performance of our benchmark methods, we use the mean squared error (MSE) on a held-out test dataset of state-action trajectories. This evaluation is conducted over ten runs, each initialized with different random seeds. We report the average MSE from these runs along with their 95% confidence intervals. Further details can be found in Appendix C. 
main results: 
 We conducted a comprehensive evaluation of our benchmark methods across all datasets, as tabulated in Table 2 . We observe that D3 can discover well-fitting dynamical system models, achieving low mean squared error in test predictions on the held-out test dataset of individual trajectories. Crucially, it can discover concise closed-form equation white-box models, of 5-15 parameters that can outperform some of the standard white-box and black-box modeling methods, indicating it is discovering well-fitting underlying equations for the respective systems of interest. , 2013] . We applied D3 to a public dataset of 33 patient trajectories [Janssen et al., 2022] , discovering a novel PK model for warfarin that outperforms existing literature (Table 3 ). We compare the test MSE of the discovered model against existing models and provide pharmacologist feedback. Full model details and results from running D3-hybrid on the same dataset are in Appendix G. Existing Wafarin PK Model. The standard pharmacokinetic (PK) warfarin model from the literature [Lv et al., 2017 , Hamberg, 2013] is the following: dC dt = k a • D -k e • C, k e = k eEQUATIONwhere the additional parameters k decay represent the natural decay rate of warfarin concentration and K m is the Michaelis constant indicating the warfarin concentration at which the metabolism rate is half its maximum. The parameters k ds , k as , and k da are the interaction terms for dosage-sex, age-sex, and dosage-age, respectively, and A and S are the sample population means for the covariates of age and sex. This new PK model includes the original parameters along with additional interaction terms and transformations, enhancing the model's complexity and accuracy. This section provides an in-depth analysis of D3's effectiveness related to its benchmark counterparts.Can D3 perform feature acquisition and leverage the LLM prior information to perform this better? To explore the feature acquisition performance of D3, we showcase this component working in full, by comparing it against the baseline of a random feature acquirer policy, a null policy that just uses the existing features and our approach of D3 that leverages the LLM to quantify the value of information for features, where no training data exists for them in the dataset. We observe, as in Figure 2 , that D3 converges the fastest, and achieves the overall highest performance, whereas the other feature acquisition methods fall short. Experimental details are in Appendix G.2.Can D3 evolve its modular model to fit the system best? We analyze this from an empirical point of view to determine if D3 can correctly evolve the generated model and reduce its prediction error over subsequent generations. We observe that D3 can indeed understand, reason, and iterate the generated code representation of the model to incorporate a better fitting model, as observed in Figure 3 . 
conclusion: 
 Summary. In this paper, we proposed Data-Driven Discovery (D3) framework, to iteratively discover and refine interpretable models of pharmacological dynamics, where it has the ability to acquire and integrate new features, validate, and compare pharmacological dynamical systems models, thereby uncovering new insights into pharmacokinetic and physiological processes. Specifically, applied to a real Warfarin PK dataset we were able to discover a better performing Warfarin PK model and provide new insights into what an optimal PK model for Warfarin may need to possess, such as additional Michaelis components.Limitations & Future work. There exist limitations to the current approach. First, the LLM discovery framework is an initial framework, and the utility of the LLM to generate better models could be enhanced, tree-based generation strategies such as the tree of thought [Yao et al., 2024] , or graph of thought [Besta et al., 2024] . Specifically, it relies on a capable enough LLM, that can use tools, where we specifically used GPT4 as the underlying LLM. Moreover, we make the assumption that when D3 acquires a new feature it acquires that feature for all the individuals within the existing dataset, we leave for future work to consider applications where this assumption may not always be true. Furthermore, we provide to the LLM a system description from a user, future works could explore avoiding the user to provide this, and instead leveraging a form of automatic retrieval augmented generation [Lewis et al., 2020] .Broader Impact and Ethical Considerations. Principally D3 aims to discover interpretable models of pharmacological dynamics. However, the final discovered models should always be checked by appropriate human experts and validated in additional held-out datasets before any clinical use. Furthermore, D3 is a tool, which could be misused by a malicious user with unethical system descriptions as input to discover a potentially biased model. Moreover, D3 leverages LLMs and their feedback as an integral component, however, LLMs are prone to hallucinations, thereby motivating any LLM-generated outputs that are shown to the user should have a content filter applied to them. Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?Answer: [Yes] Justification: We accompany results with error bars in the main results, as shown in Section 6. We explain how these are computed, and that they performed over random seed runs, which is provided in the main experimental results section, Section 6, and expanded in the evaluation metrics section, Appendix C.Guidelines:• The answer NA means that the paper does not include experiments.• The authors should answer "Yes" if the results are accompanied by error bars, confidence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper. • The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions). • The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.) • The assumptions made should be given (e.g., Normally distributed errors).• It should be clear whether the error bar is the standard deviation or the standard error of the mean. • It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have a 96% CI, if the hypothesis of Normality of errors is not verified. • For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative error rates). • If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text. 
experiments compute resources: 
 Question: For each experiment, does the paper provide sufficient information on the computer resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?Answer: [Yes] Justification: We provide full details of the compute used in Appendix C, and further detail the LLM used in Appendix F and for the respective baselines in Appendix E.Guidelines:• The answer NA means that the paper does not include experiments.• The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage. • The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute. • The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make it into the paper). 
code of ethics: 
 Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https://neurips.cc/public/EthicsGuidelines?Answer: [Yes] Justification: Yes, and we provide an explicit ethical consideration discussion section in Section 7.Guidelines:• The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.• If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics. • The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction). 
broader impacts: 
 Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?Answer: [Yes] Justification: Yes, and we provide an explicit broader impact discussion section in Section 7.Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations. • The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?Answer: [NA] Justification: It is not applicable, however, we attempt to cover such discussions in the "Broader Impact and Ethical Considerations" section in Section 7.Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters. • Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.• Including this information in the supplemental material is fine, but if the main contribution of the paper involves human subjects, then as much detail as possible should be included in the main paper. • We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution. • For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review. 
 
 Appendix: We expand on the related work contained within the main paper.ODE learning methods. Symbolic regression methods [Koza, 1994 , Kacprzyk et al., 2024 , Holt et al., 2023b] such as SINDy and D-CODE [Brunton et al., 2016 , Qian et al., 2022] can automatically discover closed form ODEs f ∈ M C through genetic algorithms. In principle, they can incorporate features a n , achieving refinement R2, but they struggle with many variables (e.g., more than 20) due to computational complexity. On the other hand, Neural ODE leverages neural networks f ∈ M N to capture a large number of variables [Chen et al., 2018 , Dupont et al., 2019 , Zaytar and El Amrani, 2016 , Devlin et al., 2018 , Sehovac and Grolinger, 2020 , Holt et al., 2022 , 2023a , 2024a] . However, Neural ODEs involve many free parameters and cannot automatically reduce the level of refinement when there is insufficient training data, leading to overfitting. In comparison, D3 can capture f ∈ M C ∪ M N and R1 -R3 levels of refinement; it can further adapt the level of refinement based on data availability. Furthermore, in both Symbolic regression and Neural ODE, human experts need to manually incorporate the problem context c through the specification of many hyperparameters. Finally, neither of these methods assumes fixed training data and does not address the Data Acquisition Challenges.AI for automated modeling. AI Feynman [Udrescu and Tegmark, 2020] improves upon the symbolic regression methods by leveraging a specific type of information c, i.e. the physical unit of each variable (e.g., meter). In applications in physics, knowing the units provides additional constraints on the variable's relationship, thus helping narrow down the search space. However, its utility outside physical sciences has been limited. More recently, Eureqa [Ma et al., 2023] applies LLMs to modeling the reward function of reinforcement learning algorithms. Eureqa is similar to D3 in its ability to integrate unstructured context c and represent the model as code. However, Eureqa cannot automatically adjust the level of refinement or acquire new data and it does not easily apply to learning temporal dynamics, e.g. ODEs. Furthermore, other LLM code generation [Holt et al., 2024b,c] approaches cannot acquire new data.Data Acquisition. Methods have been developed to address specific data acquisition challenges. Active Feature Acquisition (AFA) attempts to measure additional feature variables h on existing training samples to improve predictive performance [Ma et al., 2018 , Gong et al., 2019] . On the other hand, Active Learning [Sebastiani and Wynn, 2000 , Settles, 2009 , Sener and Savarese, 2017 , Imberg et al., 2020] and Active Testing [Lowell et al., 2018 , Kossen et al., 2021] attempt to acquire new samples to increase the size of training or evaluation data N while keeping the features constant.Neither field has proposed a method to address the zero-shot setting, where the data acquisition decision must be made before observing any target data (Equation 7). Moreover, most existing methods in Active Learning focus on acquiring labels for supervised learning given unlabelled samples, which is not immediately applicable to temporal modeling tasks.In the following, we present six clinically relevant datasets, each of which is either sourced from real-world data or generated from a high-fidelity simulator developed by domain experts.Three of our environments are derived from a state-of-the-art biomedical Pharmacokinetic-Pharmacodynamic (PKPD) model of lung cancer tumor growth, used to simulate the combined effects of chemotherapy and radiotherapy in lung cancer [Geng et al., 2017] . This model has been extensively used in other works [Seedat et al., 2022 , Bica et al., 2020 , Melnychuk et al., 2022] . We use this bio-mathematical lung cancer model to create three variations: no treatments (Lung Cancer), chemotherapy only (Lung Cancer (with Chemo.)), and both chemotherapy and radiotherapy (Lung Cancer (with Chemo. & Radio.)). For each model, we sample a respective dataset. Below, we detail the general case of Lung Cancer (with Chemo. & Radio.), which originates from the general Cancer PKPD Model, and then describe the variations.Cancer PKPD Model. This model simulates the combined effects of chemotherapy and radiotherapy on lung cancer tumor growth [Geng et al., 2017] , as shown in Equation ( 11). The model predicts tumor volume x(t) over time t in days post-diagnosis. The model includes two binary treatments: (1) radiotherapy u r t and ( 2) chemotherapyEQUATIONThe parameters K, ρ, β c , α r , β r for each simulated patient are specified in Geng et al. [2017] and summarized in Table 4 . Additionally, the chemotherapy drug concentration C(t) follows an exponential EQUATIONThe chemotherapy binary action represents increasing the C(t) concentration by 5.0mg/m 3 of Vinblastine given at time t. The radiotherapy concentration d(t) represents 2.0 Gy fractions of radiotherapy given at timestep t, where Gy is the Gray ionizing radiation dose.Time-dependent confounding. We introduce time-varying confounding by modeling chemotherapy and radiotherapy assignment as Bernoulli random variables. The probabilities p c and p r depend on tumor diameter as follows:EQUATIONwhere D max = 13cm is the maximum tumor diameter, δ c = δ r = D max /2, and D(t) is the average tumor diameter. The parameters γ c and γ r control the extent of time-varying confounding, withγ c = γ r = 2.Sampling datasets. Using the Cancer PKPD model, we sample N = 1, 000 trajectories, corresponding to N = 10, 000 patients. Initial tumor volumes are sampled from a uniform distribution x(0) ∼ U(0, 1149), and patient trajectories are forward simulated for 60 days using the Cancer PKPD Equation ( 11) and the action policy of Equation ( 13), employing a Euler stepwise solver. This process generates one dataset sample. We repeat this with independent random seeds to create D train , D val , D test . For each benchmark method run with a random seed, we re-sample the datasets. Each variation includes either chemotherapy, both chemotherapy and radiotherapy, or neither. We provide further details of this dataset's system description and variable descriptions in Appendix F.5.We utilize the sophisticated epidemic agent-based simulator COVASIM [Kerr et al., 2021] to model COVID-19 epidemics. This advanced simulator is capable of simulating both non-pharmaceutical interventions (such as lockdowns, social distancing, and school closures) and pharmaceutical interventions (such as vaccinations). In this agent-based simulator, each agent represents an individual within the population, and can be in one of several states: susceptible to COVID-19, exposed, infectious, or recovered (including deaths).We use COVASIM with its default parameter settings provided by the open-source implementation 5 . COVASIM simulates a population of individuals, and to ensure accuracy, we simulate 24 countries, collecting trajectories for each. For each simulation, we use a population size of 1, 000, 000 individuals, simulating each individual separately (disabling simulation rescaling). The simulation begins with a random number of individuals initially infected with COVID-19, I(0) = U(10, 000, 100, 000), and proceeds for 60 days.We repeat this process with independent random seeds to generate D train , D val , D test . For each benchmark method run with a different random seed, we re-sample the datasets. Detailed descriptions of this dataset's system and variables are provided in the prompt template outlined in Appendix F.5.This subsection describes an ecological model of a microcosm consisting of algae, flagellate, and rotifer populations, replicating an experimental three-species prey-predator system [Hiltunen et al., 2013] . We use the dataset provided by [Bonnaffé and Coulson, 2023] 6 . The dataset contains a single trajectory of 102 time steps. We split the data into training, validation, and test sets with proportions of 70%, 15%, and 15%, respectively, ensuring that the splits maintain the chronological order to preserve temporal causality.Detailed descriptions of this dataset's system and variables are provided in the prompt template outlined in Appendix F.5.Here we describe the real PK Wafarin dataset from a clinical trial that is publicly available [Janssen et al., 2022] . We use the NOMEN dataset available at https://github.com/Janssena/ SI-AIEP-paper. This publicly available dataset consists of 32 patients who received warfarin in a clinical trial, which was originally designed to determine how to predict drug concentrations in patients. The dataset includes a total of 251 warfarin concentration measurements, with a median of six measurements per patient. Each patient received a single dose of warfarin at t = 0, and measurements were taken at t ∈ {0.25, 0.5, 1.0, 2.0, 4.0, 6.0, 12.0, 24.0, 48.0, 72.0, 96.0, 120.0}. The available covariates in the dataset are patient weight, age, and sex. We follow the dataset's original pre-processing code, which is provided at https://github.com/Janssena/SI-AIEP-paper. We split the data into training, validation, and test sets with proportions of 70%, 15%, and 15%, respectively, ensuring that the splits maintain the chronological order to preserve temporal causality. This dataset is released under a GPL-3.0 license.Detailed descriptions of this dataset's system and variables are provided in the prompt template outlined in Appendix F.5.We utilize mean squared error (MSE) to assess the benchmark methods on a separate test dataset comprising individual trajectories, denoted as D test . This assessment is based on the loss defined in Equation ( 14) and reported as MSE. The metrics are averaged over ten runs with different random seeds, and we provide these averages along with their 95% confidence intervals. Specifically we provide the mean confidence interval7 . For each random seed, we independently generate new training, validation, and test datasets when a simulator is available. Specifically, for each simulation, we ensure that the validation and test datasets contain the same number of trajectories as the training set. Each baseline model is trained on the training dataset, with early stopping applied using the validation dataset when supported by the method. Subsequently, we evaluate the performance of each baseline on the test dataset. This entire procedure is repeated for each random seed. All experiments and training were conducted using a single Intel Core i9-12900K CPU @ 3.20GHz, 64GB RAM, and an Nvidia RTX3090 GPU with 24GB of memory.We evaluate the optimization loss using mean squared error (MSE) on a dataset D, and also consider a component-wise MSE loss for more detailed analysis.MSE Loss. Specifically, we optimize the following mean squared error objective:L(θ, D) = 1 N × T N n=1 T n ti=0 ∥f θ (x n (t i ), a n (t i ), t i )∆t -y n (t i )∥ 2 (14)where N × T represents the total number of state-action pairs in the dataset. The goal is to find the parameters θ * that minimize this loss, i.e., θ * = arg min θ L(θ, D train ). We optimize θ using stochastic gradient descent with the Adam optimizer [Kingma and Ba, 2014] , although other optimization algorithms, such as black box optimizers, could also be employed.MSE Loss per component. To gather detailed quantitative statistics on the performance of the trained system model, we collect the validation loss per component. Using (j) to denote the predictions for the j th component, we define:ω j (θ * , D val ) = 1 N val × T Nval n=1 T n val ti=0 (f θ (x n (t i ), a n (t i ), t i )∆t -y n,j (t i )) 2 (15)We aggregate these scalar validation losses per component into a vector ω = [ω 1 , ω 2 , . . . , ω m ], and compute the mean validation loss as v = 1 m m j=1 ω j (θ * , D val ).We sought feedback from three expert pharmacologists on the discovered model. They have validated and commented on the discovered model and the impact of such a tool as D3.Prof. Eoin Mckinney, Clinician. "This model is significant, as consortiums are dedicated to improving Warfarin modeling [Consortium, 2009] . The updated PK model (Equation ( 10)) adds complexity with interaction terms and introduces two new variables: the Michaelis constant and the natural decay rate of warfarin concentration. The Michaelis constant, which varies between individuals or over time due to factors like temperature and pH, is presumably estimated from observed data. The natural decay rate is likely derived from prior reports but also varies between individuals. The model has proposed minimal-cost new features (a constant and interaction terms), but it's worth noting the genotype of certain warfarin-metabolizing enzymes (CYP2C9 and VKORC1) dictates clearance rates and levels, though the model did not recommend measuring these, possibly due to the associated costs. The principal clinical relevance of the improved model is better prediction of individual warfarin doses. Numerous algorithms attempt this, supported by a consortium focused on dose-prediction models. The data used, including public datasets, supports demonstrating the clinical utility of the improved PK model, highlighting the relevance of this problem despite newer alternatives to warfarin."Jean-Baptiste Woillard, Pharmacologist. "The model is promising and pharmacokinetically plausible. I have data for mycophenolic acid, which has a challenging, complex absorption profile that would be interesting to test. The challenge extends to PK/PD, linking exposure and effect (e.g., tumor growth). Monolix, based on the SAEM algorithm, includes relevant datasets such as the PK/PD Warfarin dataset and a PK dataset for tacrolimus in heart transplants, which could serve as valuable benchmarks for comparison. Additionally, TMDD models, where clearance changes with tumor size, are of particular interest. It would be beneficial to see performance comparisons with these datasets using your approach." Richard Peck, Clinical Pharmacologist. "This model is reasonable and potentially superior. It represents a significant advance in clinical pharmacology by automatically identifying robust PK models. While I am not a modeler or an expert in warfarin modeling, the discovered model appears suitable for your purpose. The most interesting aspect of this work is that D3 identified an alternative model that could be superior. While improving PK models for warfarin might not have much clinical use since warfarin dosing is monitored by INR rather than concentration, the ability to automatically identify robust PK models is a tremendous advance. Publishing about DI3 in clinical pharmacology literature would be significant. I haven't seen anything like this before. However, to show superiority, consider studying a drug with a more challenging PK. For such a publication, involving an expert pharmacometrician would be beneficial."To evaluate the performance of D3 and establish its competitive performance, we conduct comparisons with leading modeling methods for ODEs. Specifically, we benchmark against advanced black-box models that have many parameters, such as neural ODEs with action inputs, known as DyNODE [Chen et al., 2018 , Alvarez et al., 2020] . Also in this class are recurrent neural networks (RNN) and a state-of-the-art transformer model (Transformer). In addition, we include white-box transparent dynamical systems models identified through equation discovery techniques, such as Sparse Identification of Nonlinear Dynamics (SINDy) [Brunton et al., 2016] . Moreover, D3 supports two discovery modes, of which we compare against both, discovering only white-box models (D3-white-box) and discovering hybrid models (D3-hybrid). We also perform ablations of D3, of zero-shot generated model from D3 as (ZeroShot) and the same model with optimized parameters (ZeroOptim).DyNODE DyNODE is a neural network-based model that integrates control inputs into the neural ordinary differential equation (ODE) framework [Chen et al., 2018] , as outlined by Alvarez et al. [2020] . Our implementation of DyNODE features a 3-layer Multilayer Perceptron (MLP) with hidden layers consisting of 128 units each, utilizing tanh activation functions. We initialize weights using the Xavier method [Kumar, 2017] . For consistency and competitiveness, we apply the same objective, optimizer, and hyperparameters as used in D3. Specifically, we utilize the Adam optimizer [Kingma and Ba, 2014] with a learning rate of 0.01, a batch size of 1,000, and early stopping with a patience of 20. The model is trained for 2,000 epochs to ensure convergence.RNN Recurrent Neural Networks (RNNs) [Graves et al., 2007] serve as a standard benchmark for autoregressive time series next-step prediction. In our implementation, the input dataset is normalized according to the training dataset. The model consists of a Gated Recurrent Unit (GRU) RNN that maps the state-action dimension to a hidden dimension of 250 across two layers. This hidden representation is then passed through a linear layer to convert it back to the state dimension, enabling next-step prediction. For a fair comparison, we employ the same objective, optimizer, and hyperparameters used in D3. Specifically, we utilize the Adam optimizer [Kingma and Ba, 2014] with a learning rate of 0.01, a batch size of 1,000, and early stopping with patience of 20 epochs. The model is trained for 2,000 epochs to ensure convergence.The Causal Transformer represents a cutting-edge model designed for estimating counterfactual outcomes, as detailed by [Melnychuk et al., 2022] . Given its intricate structure, which includes three distinct transformer networks for processing covariates, past treatments, and past outcomes, we opted for a streamlined approach more suitable for our datasets and task domains. Specifically, we employed a single transformer to model past outcomes. This implementation utilizes a standard transformer encoder, where the input data is normalized to the training dataset. The state-action input dimensions are encoded into a 250-dimensional embedding vector via a linear layer, followed by the application of a standard positional encoder [Melnychuk et al., 2022] . This data is then processed through a transformer encoder layer with a head size of 10 and a dropout rate of 0.1. The output from this layer is passed through a linear layer to reconstruct the next state, matching the state dimension size. Training this model involves the AdamW optimizer [Kingma and Ba, 2014] with a learning rate of 0.00005, complemented by a step learning rate scheduler with a step size of 1.0 and gamma of 0.95. Gradient clipping is applied at 0.7, and the model is trained with a batch size of 1,000. Early stopping is used with patience of 20 epochs, and training is conducted for up to 2,000 epochs to ensure convergence.SINDy Sparse Identification of Nonlinear Dynamics (SINDy) [Brunton et al., 2016] is a data-driven methodology designed to uncover the governing equations of a dynamical system directly from time-series data, resulting in a transparent, closed-form mathematical model. The SINDy algorithm operates by performing sparse regression iteratively on a library of candidate functions to find the most parsimonious and accurate representation of the system's dynamics.In our approach, we employ a second-order polynomial library, L = {1, x 0 , x 1 , x 0 x 1 }, as the feature library. To calculate the time derivatives from the input time-series data, we use finite difference approximations of first order. We maintain the alpha parameter at 0.5 for all experiments, with a sparsity threshold of 0.02, except for the COVID-19 dataset where it is set to 1 × 10 -5 . Data-Driven Discovery (D3) We defer to Appendix F for implementation-specific details. Specifically, ZeroShot and ZeroOptim are ablations of D3, where these use the exact same setup, hyperparameters, and prompts. First, ZeroShot generates one model without optimizing it's parameters θ to the training dataset split, thereby assessing the model's output loss directly from the LLM. Second, ZeroOptim re-uses the same setup as ZeroShot, however now optimizes the parameters to the data, using the training setup described in Appendix F.2.The data-driven discovery method follows the framework as outlined in Section 3. We present pseudocode in Appendix F.1, how the code-generated models f θ are trained in Appendix F.2, prompt templates in Appendix F.4, system description prompts in Appendix F.5 for each dataset. Specifically, we find a top-K, where K = 16 is sufficient. Additionally, we use the LLM of GPT4-1106-Preview, with a temperature of 0.7. f i = G(c, f i-1 , D i , r i-1 , s i-1 ) ▷ Modeling Agent G generates model f i 6: θ = arg min θ L(f i (θ), D i,train ) ▷ Fit the model f i 7:Compute validation loss L(f i ( θ), D i,val )8:s i ← s i ∪ {(f i , θ, L(f i ( θ), D val ))} ▷ Add model to top-K models 9:r ∼ E(c, f i , D i ) ▷ Generate self-reflection and feedback from Evaluation Agent E 10:if ∃ features possible to acquire then 11: Update dataset: D i,train , D i,val ← include h i 13: Return: Best model from s with the lowest validation lossh i ← A(c, f i , D i ,Upon generation by the Modelling LLM Agent, the model f is output as code, specifically a PyTorch [Paszke et al., 2017] neural network module. This code is then executed and the module is trained on the provided training dataset. The Modelling LLM Agent observes a code skeleton within the system description context c, examples of which are given in Appendix F.5. The code skeleton must be a 'torch.nn.Module' named 'StateDifferential', with initialized parameters and a forward function that computes the state differential using the state and action as input variables. The LLM completes, but does not alter, the skeleton, simplifying text processing, module execution, and model training.The model is trained using the standard MSE loss function Equation ( 14), optimized with the Adam optimizer [Kingma and Ba, 2014] , using a learning rate of 0.01, a batch size of 1,000, and early stopping with a patience of 20 epochs. Training proceeds for up to 2,000 epochs to ensure convergence and fair comparison.After training, the validation MSE and per-component MSE (validation loss per state output dimension) are computed Equation ( 15). The trained model is appended back into s i with a string representation including initialized parameter values. Providing these optimized parameters aids the LLM in suggesting effective initial values in subsequent generations, which are then refined through further optimization.Computational Efficiency. The computational efficiency of D3 overall arises from two sources: Using a Large Language Model (LLM) to generate the completions of the respective agents, and fitting the model, we now discuss each of these separately. First, using a capable enough LLM such as GPT4, for performing the completions for the agents can be computationally demanding, however, this computer costs with the computational cost per token of the underlying LLM that is used, and scales in the aspect of data with the number of tokens required to process, which we could envisage arises from having the LLM utilize a dataset with many features, each with their own textual description scaling the number of tokens in the input that the LLM has to process. One mitigation for this is to use LLMs that can handle larger context window sizes, and or be more selective about the input tokens fed into the LLM, such as restricting the number of in-context examples in the top-k examples, by reducing k, a hyperparameter. Second, we represent the models with parameters as a PyTorch module and use Pytorch to optimize the parameters, using a standard ML pipeline to train neural networks using stochastic gradient descent. Such approaches should scale with the number of input features, and parameters, however as with training any large eventual parameter model will scale as well as any other neural network-based approach in the same data pipeline. We note that more complicated distributed neural network training paradigms exist, however emphasize the focus of this work was on an initial framework, and leave the specific implementation for scale for future work.Our method is initiated by an expert who provides a structured prompt that provides the system {system description} and {feature description}.• The system description provides a holistic description of the system, specifying the state variables x(t) and actions a(t) that are used to model the system. The variables are described semantically.• The feature description define features, names, units, and their ranges.In addition to these components, the structured prompt includes a skeleton code program {skeleton program} to instruct the LLM to synthesize executable code in a pre-determined format. We refer to this prompt as the context c.For all prompt templates used, please see the code https://github.com/samholt/ DataDrivenDiscovery.By following our proposed system requirements format Appendix F.3 we constructed prompts for each of the datasets that we evaluated against, which are provided in the code, https://github. com/samholt/DataDrivenDiscovery.Hybrid PKPD Warfarin Model with Neural Network Integration. This model incorporates a simplified pharmacokinetic-pharmacodynamic (PKPD) structure with additional neural network complexity, achieving a validation loss of 0.276. The model is defined as:EQUATION) where F represents the bioavailability of the drug, V denotes the volume of distribution, k c is the clearance rate of warfarin, β a is the slope of the age effect on drug clearance, β s is the offset for the sex effect on drug clearance, MLP(C, D, A, S) is a multi-layer perceptron output, representing complex relationships in the data captured by the neural network, and R is a regularization term to prevent overfitting, with λ being the regularization weight. The addition of a neural network allows the model to learn nonlinear relationships directly from data, providing improved predictive performance.Existing Wafarin PK Model. Test MSE Loss: 0.6461. # S i m p l i f i e d MLP a r c h i t e c t u r e s e l f . mlp = nn . S e q u e n t i a l ( nn . L i n e a r ( 4 , 1 0 ) , nn . LeakyReLU ( ) , nn . L i n e a r ( 1 0 , 1 ) ) # L1 R e g u l a r i z a t i o n s e l f . r e g u l a r i z a t i o n _ w e i g h t = nn . P a r a m e t e r ( t o r c h . t e n s o r ( 0 . 0 1 ) . a b s ( ) ) d e f f o r w a r d ( s e l f , w a r f a r i n _ c o n c e n t r a t i o n : t o r c h . T e n s o r , w a r f a r i n _ d o s a g e : t o r c h . T e n s o r , p a t i e n t _ a g e : t o r c h . T e n s o r , p a t i e n t _ s e x : t o r c h . T e n s o r ) -> T u p l e [ t o r c h . T e n s o r ] : # C a l c u l a t e t h e c h a n g e i n w a r f a r i n c o n c e n t r a t i o n o v e r t i m e a b s o r p t i o n = s e l f . b i o a v a i l a b i l i t y * w a r f a r i n _ d o s a g e / s e l f . v o l u m e _ o f _ d i s t r i b u t i o n e l i m i n a t i o n = s e l f . c l e a r a n c e _ r a t e * w a r f a r i n _ c o n c e n t r a t i o n # L i n e a r e f f e c t o f a g e and s e x a g e _ e f f e c t = s e l f . a g e _ e f f e c t _ s l o p e * ( p a t i e n t _ a g e -4 0 ) s e x _ e f f e c t = s e l f . s e x _ e f f e c t _ o f f s e t * ( p a t i e n t _ s e x -0 . 5 ) a d j u s t e d _ c l e a r a n c e = e l i m i n a t i o n * ( 1 + a g e _ e f f e c t ) * ( 1 + s e x _ e f f e c t ) # MLP f o r c a p t u r i n g complex r e l a t i o n s h i p s m l p _ i n p u t = t o r c h . s t a c k ( ( w a r f a r i n _ c o n c e n t r a t i o n , w a r f a r i n _ d o s a g e , p a t i e n t _ a g e , p a t i e n t _ s e x ) , dim = 1 ) m l p _ o u t p u t = s e l f . mlp ( m l p _ i n p u t ) . s q u e e z e ( ) # R e g u l a r i z a t i o n t e r m r e g u l a r i z a t i o n = s e l f . r e g u l a r i z a t i o n _ w e i g h t * ( s e l f . c l e a r a n c e _ r a t e . a b s ( ) + s e l f . b i o a v a i l a b i l i t y . a b s ( ) + s e l f . v o l u m e _ o f _ d i s t r i b u t i o n . a b s ( ) ) d _ w a r f a r i n _ c o n c e n t r a t i o n _ _ d t = a b s o r p t i o n -a d j u s t e d _ c l e a r a n c e + m l p _ o u t p u t -r e g u l a r i z a t i o n r e t u r n ( d _ w a r f a r i n _ c o n c e n t r a t i o n _ _ d t , ) o p t i m i z e d _ p a r a m e t e r s = { ' c l e a r a n c e _ r a t e ' : 0 . 0 8 0 1 6 5 3 9 3 6 5 0 5 3 1 7 7 , ' b i o a v a i l a b i l i t y ' : 3 . 1 7 0 1 4 8 1 3 4 2 3 1 5 6 7 4 , ' v o l u m e _ o f _ d i s t r i b u t i o n ' : 3 4 . 7 4 9 3 2 8 6 1 3 2 8 1 2 5 , ' a g e _ e f f e c t _ s l o p e ' : -0 . 0 7 7 9 1 9 5 2 0 4 3 7 7 1 7 4 4 , ' s e x _ e f f e c t _ o f f s e t ' : -0 . 1 6 8 0 9 0 4 3 2 8 8 2 3 0 8 9 6 , ' r e g u l a r i z a t i o n _ w e i g h t ' : -0 . 0 6 0 7 5 4 0 9 0 5 4 7 5 6 1 6 4 6 } w h i t e and b l a c k box model w i t h s i m p l i f i e d l i n e a r a g e and s e x e f f e c t s , a s i m p l i f i e d MLP a r c h i t e c t u r e f o r c a p t u r i n g complex r e l a t i o n s h i p s , and L1 r e g u l a r i z a t i o n t o p r e v e n t o v e r f i t t i n g .For this experiment, we adapted the Lung Cancer (with Chemo. & Radio.) dataset simulator, to form a new variation of this dataset, where D3 starts with two features, which are the main state x(t) features that it intends to model, and can optionally acquire new features a(t) from a set of possible 22 features, which are listed below. Of importance, two features directly affect the underlying behaviour of this dynamical system, these being radiotherapy_dosage and chemotherapy_dosage, and the rest of the features we model as random white noise features, hence have no effect on the underling system, however D3 is unaware of this, and must use the feature description, from the features name to determine relative value of improved performance that it can attain by acquiring that feature, and propose new features to acquire along with their expected improvement.Specifically, in this setup, we allow the feature acquistion agent to select from the following set of features at each iteration:[ b l o o d _ t y p e , s m o k i n g _ s t a t u s , p a t i e n t _ a g e , p a t i e n t _ w e i g h t , p a t i e n t _ h e i g h t , sex , We also make the assumption that each feature that can be acquired has equal cost l of acquisition, and that the agent has an unlimited feature acquisition budget, however these settings can be readily changed depending on the specific application.To provide comparable baselines, we compare D3 against ablations of itself, which is a random feature acquirer Feature Acquisition Agent, that selects a feature randomly for acquisition, and a null, or no Feature Acquisition Agent, that only observes at all iterations the target states x(t) to model. We observe experimentally that D3 that can leverage the LLM to quantify the value of information for features, where no training data exists for them in the dataset. We observe, as in Figure 2 , that D3 converges the fastest, and achieves the overall highest performance, whereas the other feature acquisition methods fall short.We provide part of the logs from one of the random seed runs, including part of it due to the logs being extensive. The t i m e u n i t s i s i n d a y s . 2 5 8 5 7 0 4 2 5 2 1 2 3 8 3 2 7 , ' b e t a ' : 0 . 0 2 5 8 6 2 5 1 8 6 9 7 9 7 7 0 6 6 , 'gamma ' : 0 . 0 3 7 3 0 2 4 3 8 1 6 9 7 1 7 7 7 0 4 , 7 0 In the main paper we presented the best Warfarin discovered white-box model, however after running D3-white-box across all the other datasets, we include a few of the white-box models that it discovered and generated for each other respective dataset used that we compared against. Of particular interest, is that these white-box models fit the dataset well, and can even outperform some of the other neural network based baselines-we highlight that with the right functional form of the white-box model, we can achieve better performance than a largely over-parameterized neural network model, that has signficiantly more parameters (such as an RNN with 18,000 parameters), than a white-box model with 5-15 parameters.Val L o s s : 0 . 2 6 5 7 3 6 8 5 4 6 3 5 1 7 9 0 4 , ' b e t a ' : 1 0 0 0 . 9 4 5 9 8 3 8 8 6 7 1 8 8 , 'gamma ' : -0 . 0 0 7 2 1 9 3 5 6 0 1 9 0 4 9 8 8 3 p a r a m e t e r s = { ' a l p h a ' : 0 . 0 4 9 9 9 9 9 9 3 2 9 4 4 7 7 4 6 , ' b e t a ' : 0 . 2 3 9 9 9 9 9 9 4 6 3 5 5 8 1 9 7 , 'gamma ' : 0 . 2 3 9 9 9 9 9 2 0 1 2 9 7 7 6 , ' d e l t a ' : 0 . 1 3 9 9 9 9 9 8 5 6 9 4 8 8 5 2 5 , ' e t a ' : 0 . 0 5 0 0 0 0 3 0 9 9 4 4 1 5 2 8 3 , ' t h e t a ' : -0 . 1 2 9 9 9 9 2 7 9 9 7 5 8 9 1 1 , ' kappa ' : 2 . 2 0 9 9 9 9 7 9 9 7 2 8 3 9 3 6 , ' sigma ' : 1 . 3 9 0 0 0 0 7 0 0 9 5 0 6 2 2 6 , ' rho ' : 0 . 6 9 0 0 0 0 2 9 5 6 3 9 0 3 8 1 , ' omega ' : 0 . 0 1 9 9 9 9 9 9 9 5 5 2 9 6 5 1 6 4 , ' p s i ' : 0 . 1 2 9 9 9 9 9 5 0 5 2 8 1 4 4 8 4 , ' p h i ' : 0 . 0 2 0 0 0 0 0 0 3 2 7 8 2 5 5 4 6 3 , ' x i ' : 0 . 0 2 0 0 0 0 0 9 0 8 2 2 5 7 7 4 7 7 } ### NeurIPS Paper Checklist• The answer NA means that the paper does not include experiments.• If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not. • If the contribution is a dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable. • Depending on the contribution, reproducibility can be accomplished in various ways.For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution is a specific model and empirical evaluation, it may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., in the case of a large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed. • While NeurIPS does not require releasing code, the conference does require all submissions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example• The experimental setting should be presented in the core of the paper to a level of detail that is necessary to appreciate the results and make sense of them. • The full details can be provided either with the code, in appendix, or as supplemental material.Code is available at https://github.com/samholt/DataDrivenDiscovery and we provide a broader research group code base at https://github.com/vanderschaarlab/DataDrivenDiscovery.We provide full commentary in Appendix D.COVASIM is an open-source simulator, available at https://github.com/ InstituteforDiseaseModeling/covasim.The Plankton Microcosm and Hare-Lynx datasets are both open source and available at https://github. com/WillemBonnaffe/NODEBNGM.We use the following code https://stackoverflow.com/questions/15033511/ compute-a-confidence-interval-from-sample-data. 


Decision: As the Area Chair, I have carefully reviewed the submission, the reviewers' comments, and the authors' rebuttal. The paper presents a novel framework called Data-Driven Discovery (D3) that leverages Large Language Models (LLMs) to discover and refine interpretable models of pharmacological dynamics. The reviewers generally agree that the approach is innovative and shows promise in the field of pharmacokinetic modeling.

The main strengths of the paper include its well-structured presentation, clear methodology, and the framework's ability to iteratively improve models and acquire new features. The authors demonstrated the effectiveness of their approach on both simulated benchmarks and a real-world Warfarin dataset, which adds to the practical relevance of the work.

However, the reviewers raised several concerns. One reviewer pointed out the lack of comparison with recent symbolic regression methods like D-CODE. In response, the authors clarified that their focus was on proposing a method that can discover interpretable models, incorporate textual prior, and acquire new features on demand, which existing symbolic regression methods cannot do.

Another concern was the potential for knowledge leakage due to pre-trained LLMs. The authors addressed this by conducting additional experiments on five semi-synthetic datasets, demonstrating that D3-Hybrid performs well even on models the LLM has never seen before. This provides evidence to mitigate the concern of knowledge leakage.

The reviewers also noted the absence of ablation studies and case studies on different module choices. The authors clarified that they did include ablations in most experimental results, specifically with the additional baselines of ZeroShot and ZeroOptim models. They have revised the baseline method description to make these ablations more prominent.

Some concerns were raised about the scalability, computational efficiency, and potential for hallucinations in LLMs. While the authors acknowledged these limitations, they did not fully address how they plan to mitigate these risks in future work.

In conclusion, while there are some limitations and areas for improvement, the paper presents a novel and promising approach to pharmacokinetic modeling using LLMs. The authors have addressed many of the reviewers' concerns in their rebuttal and provided additional experimental results to support their claims. 

Based on the overall assessment, I recommend accepting this paper for publication at NeurIPS. The innovative use of LLMs in pharmacokinetic modeling, combined with the authors' thorough response to reviewer concerns, makes this work a valuable contribution to the field. However, I suggest that the authors incorporate their additional experiments and clarifications into the final version of the paper to address the reviewers' concerns more comprehensively.

Reviews: review: paper_summary: paper_summary: This paper develops an LLM-assisted equation discovery framework especially for pharmacokinetic process. Three agents of Modeling Agent, Evaluator Agent and Feature Acquisition Agent are built to explore, refine and iterate vast model space, including three levels of initial conditions, observed features and possible acquired features. Several experiments mainly on simulated benchmarks are conducted to evaluate the performance of the proposed framework. Overall, this paper investigated whether LLM agents can help effectively search and determine the model space for dynamics modeling and process discovery in pharmacology.
strength_weakness: strength_weakness: 1. The benchmark methods are not recent works, which might not well identify the effectiveness of the proposed framework. For example, the symbolic regression methods like D-CODE [ICLR'22] are not compared. 
2. To my understanding, the D3 model as well as other symbolic regression models discover equations from data-driven training, and then evaluate the performance based on the discovered equation. If I understand correctly, to me the experiments are not enough to prove the model effectiveness because of only five datasets with five equations. I would like to set several simulated datasets with different simulation equations for further investigation.
3. There is no ablation study on model design or case study on different module choice, which makes it hard to evaluate the robustness of the proposed framework. Also, can LLM agents provide explanations on their output results to make a more transparent solution.
questions: questions: 1. How is the memory module implemented in this paper, and what exactly is the memory $s_i$? Are there several choices for memory and how would they affect the performance?
2. How much is the cost for training a D3 framework on certain datasets, especially how much would it cost for calling GPT API? Can other open-sourced LLM be used for same tasks and how is their performance?
3. RNN and Transformer models usually perform similarly in prediction tasks. However, in Table 2, RNN and Transformer perform differently on several datasets. For example, on Lung Cancer, RNN's MSE is 1.16e6, which I am wondering if the training converges successfully.
limitations: limitations: The pharmacologists feedback statements highly depend on the personal knowledge and scope, as well as their understanding to the dataset or modeling. In my opinion, this part especially with only three experts are not worth for evaluation and should be excluded to avoid bias.
paper_summary: paper_summary: The paper proposes the Data-Driven Discovery (D3) framework, which leverages Large Language Models (LLMs) to iteratively discover and refine interpretable models of pharmacological dynamics. This approach enables the LLM to propose, acquire, and integrate new features, validate, and compare pharmacological dynamical systems models. The framework is demonstrated on a real pharmacokinetic datasets, highlighting its potential for clinical applications.
strength_weakness: strength_weakness: 1. Limited Novelty: The paper has very limited novelty in terms of machine learning and does not propose any new tools. It effectively uses existing toolboxes for a clinically relevant task. While applying an existing toolbox in the context of a new task is perfectly fine, similar methods have been published before, particularly those related to automated science labs(self driving labs).
2. Lack of Guardrails for Hallucinations: Models like GPT-4 are prone to hallucinations and can generate incorrect information. The framework does not appear to have guardrails to mitigate this risk, which is a serious concern given the clinical relevance of the task. While the authors have acknowledged this as a limitation, it remains a critical issue.
3. Prompt Dependency: The model’s performance is highly dependent on the prompts provided. It is crucial to formulate the prompts correctly for the LLM to generate reasonable responses. This dependency suggests that the prompting framework might need adjustments for each dataset, leaving users vulnerable to the LLM’s unpredictability.
4. Scalability and Computational Efficiency: The framework’s scalability and computational efficiency are not thoroughly discussed. Given the iterative nature of the model refinement process, it could be computationally intensive compared to other baselines and may not be practical for larger datasets without significant computational resources.
5. Better performance of other baselines : The performance of the D3 framework might be dependent on the data. The paper shows better performance of models like RNN and transformers on tasks such as COVID-19 and Warfarin.  The paper does not any comment on why that might be the case.
questions: questions: 1. Given that models like GPT-4 can generate incorrect information, what measures have you considered or could implement to mitigate this risk, especially for clinically relevant tasks? I think implementing a RAG like model to address hallucinations will be highly relevant in this case and make the paper better. 
2. How do you ensure consistency and reliability in the prompts across different datasets? Have you considered strategies to standardize the prompting process to minimize variability and improve predictability? Did you experiment with prompts multiple times to arrive at a reasonable framework?
3. Can you provide more information on the scalability and computational efficiency of your framework? How does the iterative model refinement process impact computational resources compared to other baselines, and what are the practical implications for larger datasets ?
4. The paper shows that models like RNN and transformers perform better on tasks such as COVID-19 and Warfarin. Can you explain why this might be the case and what limitations the D3 framework has in these scenarios? How might the performance be improved?

If the authors address some of the concerns and reevaluate the frameworks with the clinicians, I will be happy to change the score. I get a sense from the expert comments that the datasets evaluated upon might not be challenging enough.
limitations: limitations: The authors have acknowledged the limitations of their model. However, I think some of the points mentioned, such as retrieval-augmented generation, must be implemented in the framework to make the model more robust.
paper_summary: paper_summary: The paper presents the Data-Driven Discovery (D3) framework, a novel approach that iteratively discovers and refines interpretable pharmacological dynamical models using LLMs. This framework is novel and innovative in its domain, it is designed to address limitations in traditional pharmacokinetic (PK) modelling, which often relies on human expertise, and prior data in various formats. It is also constrained by existing knowledge. D3 leverages LLMs to propose and integrate new features, validate suggested models, and uncover new insights in PK. The framework demonstrates its effectiveness through experiments on various PK datasets, including a real-world Warfarin dataset, where it identified a new, well-fitting PK model outperforming existing literature.  The D#  framework utilises the following agents:
Modelling agent: Writes Python code for an AI model with the information the LLM has acquired    
Evaluator agent: Evaluates the performance of the model by the modelling agent. The metric used is MSE (Mean Squared Error) 
Feature acquisition agent: Based on the evaluation of the model, more features are added to the model to improve accuracy (reduce MSE). A feature is acquired based on an estimation of the feature’s value using existing frameworks and available information such as text-based description of the feature, summary statistics and the feature’s existing data.
strength_weakness: strength_weakness: In general, the paper is inspiring and is solid in its experimentations and validations. There are a few remarks, some of them for future consideration, while others might improve the clarity of the paper:

1. Cost function considerations: the current description of the cost function l(h_i)  might be too simplistic and not fully capture the complexity of acquiring new features. Cost functions are often difficult to quantify accurately. Adding a few sentences or a small subsection to discuss the aspects of the cost function and its potential impact and tradeoffs will provide a more comprehensive understanding. 
2. The D3 framework assumes that once a feature is selected, it is available for all individuals in the dataset. For example, certain biomarkers can be measured only in a small sample of patients. It would be interesting to discuss potential performance and limitations in light of missing data. Note: the authors briefly outlined it but elaboration on how to address it would be good given the commonality of this scenario. 
3. While authors called out the scalability and computational complexity of the framework, it would be interesting to see specific numbers around time to completion wrt the numbers of features. 
4. In line 212 it is mentioned that ‘The Evaluation Agent dynamically assesses both model performance and plausibility <…>’. Not sure if understood how plausibility is measured (couldn’t find it in the prompt in the Appendix either). Would be nice to see a bit more elaboration.
5. Table 2 has mixed formatting: mixing scientific and regular notation: like 0.000245 and 2.47e+03 are both present. Unifying the notation will be more consistent. In addition, the variability of the results in terms of MSE is concerning. Might be useful to think of a more coherent metric or to provide an explanation for such extreme variability (e.g. 5780 vs. 7.07 for Lung Cancer or 719 vs. 0.3 for Lung Cancer with Chemo ).   
6. The choice of LLM is quite significant and it is stated for the first time in Limitations on page 9. This could be brought earlier.
questions: questions: * Why was MSE used as the metric for accuracy?  And on the same topic, in Table 2, what would be an interpretation of units in MSE? 
* Do you have ideas as to how to combat the limitations mentioned such as hallucinations? 
* Do you have a view on how to address the cost/ability of the LLM to distinguish between “easier” and “harder” features?  Such as: “don’t include biomarker x as it’s difficult to acquire from patient”? 
* Have you considered the inverse bias problem when models such as GPT-4 are restricted not to use biased features but it could be informative in the pharmaceutical environment to a PK model, such as ethnicity? 
* I'm still not fully clear how the LLM acquires the unstructured context/metadata about each feature.
limitations: limitations: Outlined in the previous sections. They are outlined and very helpful but some of them spark interesting questions and would like to see some of them elaborated on a bit more like missing data.
paper_summary: paper_summary: This paper presents the D3 Data Driven Discovery framework which uses GPT4-1106-Preview in a framework to iteratively consider modifying the features used in the ODE.  The dynamical systems are evaluated on MSEs of held-out state-action trajectories.
strength_weakness: strength_weakness: The model space for the Modeling Agent could be better scoped---it appears to be whatever code the GPT-4 response is to the query.  What happens in event of failure?  Is this an automated framework?

Typically you want interpretable and performance pharmacokinetic models, so there is a tradeoff between validation performance (e.g. MSE) and simplicity (measure in some way, e.g. MDL, number of features, etc), which creates a frontier for each model class.  Could this be another way to compare the performance of these models?  An analogous question would be related to acquisition costs $l(h_i)$ as well.
questions: questions: What happens when the Modeling Agent produces an invalid model? What models are possible?
limitations: limitations: There is a limitation section, though the limitations appear more to be feature requests, rather than considerations that users would need to be aware of in using or building upon the work.
paper_summary: paper_summary: The paper  presents the data discovery framework. It uses a LLM to iteratively discover and refine interpretable models of pharmacological dynamics. The D3 framework consists of three agents, which work collaboratively to propose, acquire, and integrate new features, validate models, and uncover new insights into pharmacokinetic processes. The framework is designed to address the limitations of traditional pharmacokinetic models that are often constrained by human expertise and scalability issues. The D3 framework was validated using a real pharmacokinetic dataset of Warfarin, demonstrating its ability to uncover a new, plausible pharmacokinetic model that outperforms existing models.
strength_weakness: strength_weakness: All simulated datasets (and the ODE parameters) used in this study were publicly available before GPT-4's knowledge cutoff date. It is likely that OpenAI trained GPT-4 on relevant literature, leading to a significant knowledge leakage issue. I appreciate the author's perspective that LLMs encompass a vast amount of potentially usable knowledge. However, I do not believe that allowing LLMs to see the standard answers during training is an appropriate experimental setting. It will be better If the authors could demonstrate the feasibility of this method using real complex high-dimensional datasets whose dynamic is really unknown. After all, this is the kind of task faced in real pharmacokinetic studies. Meanwhile, the difference between this paper and many previous pharmacokinetic studies is that this paper does not conduct experiments on causal discovery. Instead, the model performance is evaluated through the MSE. It seems that this paper does not actually need to rely on simulated data and could use real high-dimensional time-series observational data for experiments.

The Table 3 did not report the performance of baseline models. I believe that comparing only the so-called "standard model" in Table 3 is inappropriate. Researchers in many other disciplines do not prioritize performance to the same extent as those in the machine learning community. They often trade off performance for model simplicity, valuing clear and concise formulas over marginal gains in performance. Thus, it is not very convincing that the authors achieved a performance advantage by using a method that resulted in a much more complex system compared to the standard model.

The baselines used in this paper are somewhat outdated. More recent studies published in the last two years, such as D-CODE and PGM (arxiv 2105.02522), could perform similar tasks as this paper. It is unclear why these were not included in the comparisons.
questions: questions: NA
limitations: limitations: NA

Comments: title: Concerns on the evaluation
Concerning the frontier
Gratitude for Revised Review and Score Increase
comment: Thank you for your prompt feedback and for updating your score in response to our clarifications on guardrails and hallucinations. We greatly appreciate your thoughtful consideration of our responses.

**Regarding MSE and Interpretability:**

We understand your concerns about the MSE differences between the D3-white box model and the transformer across datasets. While it’s true that the D3-white box model may not match the transformer in terms of raw MSE, we’d like to emphasize a few key points:

- **Interpretability’s Value in Clinical Contexts:** Even with a higher MSE, the interpretability of the D3-white box model provides crucial insights into the underlying pharmacological processes. This can lead to more informed clinical decisions, which may not be possible with black-box models, regardless of their lower MSE.

- **D3-Hybrid’s Balanced Approach:** Importantly, we recommend the D3-Hybrid model for practical use. The D3-Hybrid model combines the interpretability of the white-box approach with the performance benefits of data-driven models. It achieves a strong balance, with an average of 245 parameters—significantly fewer than the 2,558,348 parameters in a typical transformer—while still delivering competitive performance. This model retains a largely interpretable component, which is critical for clinical applications, ensuring that practitioners can trust and understand the results without sacrificing too much in terms of accuracy.

Given these considerations, we respectfully ask if you might reconsider your score once more, recognizing the unique value that the D3-Hybrid model brings. It offers a balanced, practical solution with substantial interpretability and strong performance, particularly in settings where understanding the model’s behavior is as important as its predictive accuracy.

Thank you again for your time and continued engagement.
I don't think the gains in MSE are marginal when comparing the performance of the transformer versus the D3-white box model in Table 2. The MSE of the D3-white box is roughly 5-10 times higher than that of the transformer across almost all datasets. The point was, if a model is already a poor fit, is its interpretability reliable?

I have also updated my score in response to the comments on lack of guardrails and hallucinations as the authors pointed out they did not observe any significant deviation from expert consensus
Thank you for your continued engagement and constructive feedback. We appreciate your insights, which have helped us further clarify our contributions.

**Guardrails and Hallucinations:**

We understand your concerns about the potential for LLM-generated hallucinations, particularly in a clinical setting. To address this:

1. **Validation Feedback Loop:** D3’s iterative process rigorously evaluates each model against a validation dataset, ensuring that only empirically sound models are refined and considered.

2. **Biological Plausibility:** The LLM is guided to generate biologically plausible models by incorporating prior knowledge. While we did not observe any cases of significant deviation from expert consensus, we’ve implemented a human-in-the-loop process to verify all models before finalization.

3. **Future Enhancements:** We acknowledge the potential of integrating retrieval-augmented generation (RAG) techniques to further reduce the risk of hallucinations and will discuss this as a future direction in the final paper.

**Model Interpretability vs. Performance:**

While D3’s white-box models may not always match the predictive accuracy of black-box models, they offer critical interpretability, which is essential in clinical applications. The ability to understand and trust these models is often more valuable than marginal gains in MSE. Additionally, D3’s emphasis on scalability and efficiency makes it well-suited for clinical use, where computational resources may be limited.

**Final Remarks:**

We hope this response clarifies our contributions and addresses your concerns. We kindly request you to reconsider your score in light of these points and remain open to further discussion if needed.

Thank you for your thoughtful review.
The authors mention, “While black-box models like RNNs may perform well in MSE, D3’s strength lies in offering interpretable models that provide valuable clinical insights, balancing performance with the necessity for interpretability.”

However, the D3 white-box models, which should be interpretable, show performance levels that are extremely poor compared to the black-box models. This raises concerns about whether the interpretability of such models is truly reliable. I believe a certain level of performance is necessary to ensure that the model’s interpretability is actually useful. In Table 2, the white-box model consistently performs poorly.

Another concern is that the authors claim such a model can be useful in clinical settings, but the lack of guardrails is worrying. Although the authors have attempted to mitigate the risks by incorporating a human in the loop and warning readers about potential risks, I am curious whether they observed cases where the model predicted something completely different from what human experts would agree with.
Thank you for the detailed feedback and for highlighting areas where additional clarification would be beneficial. We appreciate your focus on the balance between model performance (e.g., MSE) and model complexity (e.g., parameter count), particularly in the context of pharmacokinetic modeling. 

### Addressing Model Complexity and Performance Trade-offs

We understand the reviewer's concern about determining the frontier of model performance, especially when comparing the complexity of neural network-based models (black-box) with interpretable models (white-box). To address this, we conducted a thorough analysis of the average parameter count (running each baseline across ten random seeds) as a proxy for model complexity, and the results are summarized in the table below:

|Baseline|Parameter Count|Warfarin PK MSE $\downarrow$|
|---|---|---|
|DyNODE|33,922|0.726$\pm$0.17|
|SINDy|13|6.84 $\pm$ 1.76| 
|RNN|569,002|0.0495 $\pm$ 0.0406|
|Transformer|2,558,348|1.33 $\pm$ 0.941|
|D3-white-box|8| 19.6 $\pm$ 40.3|
|D3-hybrid|245|0.647 $\pm$ 0.167|

In our response, we chose parameter count as the metric to compare the complexity of different models for several reasons:

1. **Interpretability**: Parameter count provides a straightforward metric to quantify model complexity, which is particularly relevant when discussing trade-offs between performance and interpretability. Fewer parameters generally imply a simpler, more interpretable model, which is crucial in pharmacokinetics where model transparency is valued.

2. **Consistency**: This metric allows for a fair comparison across different types of models—ranging from highly parameterized black-box models to more concise white-box models. For instance, the D3-white-box model, despite having a significantly lower parameter count (8 parameters), still delivers competitive performance, albeit with higher variance, which reflects its simplicity and interpretability.

3. **Practical Relevance**: In pharmacokinetics, overly complex models (e.g., RNNs and Transformers) might offer better fit (lower MSE), but their high parameter count can obscure the underlying biological mechanisms, making them less useful for clinical interpretation and decision-making.

---

### Conclusion and Final Remarks

We believe that the additional analysis and clarifications provided here address the key concerns raised in your review. We hope that these points demonstrate the robustness of our approach and the careful consideration we have given to the trade-offs between model complexity and performance. We would be grateful if you could reconsider your score in light of these clarifications. Of course, we are open to further discussion should you have any additional questions or require further details.

Thank you again for your thoughtful review and for the opportunity to clarify these important points.
Thank you for the constructive feedback. Below, we address your main concerns briefly:

1. **Novelty:** While D3 builds on existing methods, it uniquely integrates LLMs for dynamic feature acquisition and iterative model refinement, which is not present in existing frameworks like automated science labs. This allows us to uncover new, interpretable models tailored for clinical applications.

2. **Hallucinations:** We mitigate hallucination risks with a human-in-the-loop process where all LLM outputs are validated by clinicians. We acknowledge the suggestion of using RAG techniques and will include a discussion in the final paper to explore this as a future direction.

3. **Prompt Dependency:** D3 addresses prompt variability through iterative refinement, ensuring prompts are well-tuned for each dataset, thus enhancing consistency and predictability across different datasets.

4. **Scalability:** Despite the iterative nature, D3 is computationally efficient due to its focus on discovering models with fewer parameters, making it practical for larger datasets while maintaining interpretability, which is critical in clinical settings.

5. **Performance of Baselines:** While black-box models like RNNs may perform well in MSE, D3’s strength lies in offering interpretable models that provide valuable clinical insights, balancing performance with the necessity for interpretability.

We hope this clarifies our contributions and addresses your concerns. We kindly request you to reconsider your score based on this summary. We are open to further discussions if needed.

Thank you for your time.
We sincerely thank the reviewer for their time and careful review of our submission. We are pleased to see that the reviewer found no ethical issues with our research and acknowledged that the study does not involve human subjects, which aligns with our stated description in the paper.

We appreciate the reviewer's careful consideration and their recommendation, and we would like to take this opportunity to reiterate the key contributions of our work:

- **Novelty and Impact**: Our paper presents the \textit{Data-Driven Discovery (D3)} framework, which is a novel approach leveraging Large Language Models (LLMs) to iteratively discover and refine interpretable models of pharmacological dynamics. We believe this contribution is significant as it has the potential to advance pharmacokinetic modeling and offer new insights into the dynamics of pharmacological systems.

- **Validation and Results**: We demonstrated the effectiveness of D3 through experiments on diverse pharmacokinetic datasets, including a real-world Warfarin dataset. The results indicate that our approach can uncover new, plausible pharmacokinetic models that outperform existing literature, particularly in the context of precision dosing in clinical applications.

- **Ethical Considerations**: We acknowledge the importance of ethical considerations in research, particularly in fields with potential clinical applications. Although our work does not involve human subjects, we remain committed to ensuring that our methods are used responsibly and that any models developed using D3 are rigorously validated before being applied in real-world clinical settings.

We hope that these clarifications and highlights further demonstrate the value and relevance of our work. We are committed to advancing this line of research and ensuring its impact is both positive and ethically sound.

---

*We hope that most of the reviewer’s concerns have been addressed. We’d be happy to engage in further discussions.*
**We thank the reviewer for their time and careful review.**

We appreciate your positive feedback regarding our work on using large language models (LLMs) to discover interpretable pharmacokinetic models. We are pleased to hear that you found our acknowledgment of the potential societal impacts sufficient and that our approach to caution readers about the risks associated with applying these models without human evaluation was appropriate.

In response to your comments, we would like to reiterate and clarify the following points:

- **Ethical Issues:** You correctly noted that our work does not involve human subjects or data, which minimizes direct ethical concerns. However, we recognize the broader ethical implications of deploying AI-discovered pharmacokinetic models in clinical settings. As we emphasized in the paper, we strongly caution against using these models without rigorous human expert evaluation. The safety and well-being of patients are paramount, and our framework is designed to assist, not replace, human judgment.

- **Potential Harm Without Human Evaluation:** We agree that models discovered by LLMs, when used without proper human oversight, could lead to harm. This is why we have explicitly stated in our paper that the models generated by our framework should always be validated by experts in the field before clinical application. We aim to develop tools that augment human expertise rather than replace it, ensuring that our models serve as valuable starting points for human-driven refinement and validation.

- **Acknowledgment of Societal Impacts:** We are glad that you found our discussion on the societal impacts to be adequate. We believe that transparency in AI research, especially in sensitive domains like healthcare, is critical. We have made a conscious effort to highlight both the benefits and risks of our approach, ensuring that readers are aware of the importance of responsible deployment.

In summary, we are grateful for your positive assessment of our work and your acknowledgment that we have sufficiently addressed the potential societal impacts. We are committed to advancing AI-driven discovery while maintaining the highest ethical standards, and your feedback has been invaluable in helping us refine our presentation.

---

*We hope that most of your concerns have been addressed. We’d be happy to engage in further discussions or provide additional clarifications if needed.*
I still have some concerns about determining the frontier of model performance when providing a small interpretable model (white-box) and a performant model (black-box), and then describing their performance respectively with "new pharmacokinetic model" and "low MSE" respectively.

Table 2 suggests D3-hybrid is similar to or marginally better than "transformer" (MSE lower in 4/6 datasets), and that D3-white-box is similar to or marginally better than SINDy/ZeroShot (though this latter comparison is very hard because we are not quantifying how sparse these white-box models are).

More compelling would be to train each model across one or more hyperparameters and plot performance curves across performance (MSE) and some form of sparsity (MDL, etc).
Dear Reviewer 2i8U,  
I am a NeurIPS 2024 Area Chair of the paper that you reviewed.

This is a reminder that authors left rebuttals for your review. 
We need your follow up answers on that. Please leave comment for any un-answered questions you had, or how you think about the author's rebuttal. 
The author-reviewer discussion is closed on Aug 13 11:59pm AoE.

Best regards,
AC
Dear Reviewer QDKe, 
I am a NeurIPS 2024 Area Chair of the paper that you reviewed.

This is a reminder that authors left rebuttals for your review. 
We need your follow up answers on that. Please leave comment for any un-answered questions you had, or how you think about the author's rebuttal. 
The author-reviewer discussion is closed on Aug 13 11:59pm AoE.

Best regards,
AC
Dear Reviewer b77T, 
I am a NeurIPS 2024 Area Chair of the paper that you reviewed.

This is a reminder that authors left rebuttals for your review. 
We need your follow up answers on that. Please leave comment for any un-answered questions you had, or how you think about the author's rebuttal. 
The author-reviewer discussion is closed on Aug 13 11:59pm AoE.

Best regards,
AC
I thank authors for their thoughtful response, I keep my ratings as-is but the improvements do strengthen the paper. Thank you!
We greatly appreciate your thoughtful review and the time you invested in evaluating our paper. Your feedback was crucial in refining our work, leading to additional experiments and improved explanations. We are thankful for your positive reassessment and the increased score. Thank you again!
Thanks for the authors' reply. I appreciate that and would like to raise my score.
Thank you very much for your thoughtful consideration and the time you have dedicated to reviewing our paper. Your feedback was instrumental in enhancing our work, through extensive new experiments and explanations, and we are grateful for your increased score. Thank you once again!
I appreciate the response from authors and I have improved my score.