Automatically Learning Hybrid Digital Twins of Dynamical Systems
Authors: Samuel Holt
Tennison Liu
Mihaela van der Schaar
Venue: NeurIPS 2024 spotlight (2024)
Keywords: Hybrid Models
Dynamical Systems
Primary area: machine_learning_for_other_sciences_and_fields

Abstract: Digital Twins (DTs) are computational models that simulate the states and temporal dynamics of real-world systems, playing a crucial role in prediction, understanding, and decision-making across diverse domains. However, existing approaches to DTs often struggle to generalize to unseen conditions in data-scarce settings, a crucial requirement for such models. To address these limitations, our work begins by establishing the essential desiderata for effective DTs. Hybrid Digital Twins (**HDTwins**) represent a promising approach to address these requirements, modeling systems using a composition of both mechanistic and neural components. This hybrid architecture simultaneously leverages (partial) domain knowledge and neural network expressiveness to enhance generalization, with its modular design facilitating improved evolvability. While existing hybrid models rely on expert-specified architectures with only parameters optimized on data, *automatically* specifying and optimizing HDTwins remains intractable due to the complex search space and the need for flexible integration of domain priors. To overcome this complexity, we propose an evolutionary algorithm (**HDTwinGen**) that employs Large Language Models (LLMs) to autonomously propose, evaluate, and optimize HDTwins. Specifically, LLMs iteratively generate novel model specifications, while offline tools are employed to optimize emitted parameters. Correspondingly, proposed models are evaluated and evolved based on targeted feedback, enabling the discovery of increasingly effective hybrid models. Our empirical results reveal that HDTwinGen produces generalizable, sample-efficient, and evolvable models, significantly advancing DTs' efficacy in real-world applications.

Automatically Learning Hybrid Digital Twins of Dynamical Systems
Samuel Holt, Tennison Liu' & Mihaela van der Schaar
DAMTP, University of Cambridge
Cambridge, UK {sih31, t1522, mv472)&cam.ar ac.aammacc... ac.uk
Abstract
Digital Twins (DTs) are computational models that simulate the states and temporal dynamics of real-world systems, playing crucial role in prediction, understanding. and decision-making across diverse domains. However, existing approaches to unseen conditions in data-scarce settings,
DTs often struggle to generalize a crucial requirement for such models. To address these limitations, our work begins by establishing the essential desiderata for effective DTs. Hybrid Digital Twins (HDTwins) represent a promising approach to address these requirements, modeling systems using a composition of both mechanistic and neural components. This hybrid architecture simultaneously leverages (partial) domain knowledge and neural network expressiveness to enhance generalization, with modular design facilitating improved evolvability. While existing hybrid models rely expert- specified architectures with only parameters optimized on data, automatically specifying and optimizing HDTwins remains intractable due to the complex search space and the need for flexible integration domain priors. To overcome this complexity, we propose evolutionary algorithm (HDTTinGen)) that employs Large Language Models (LLMs) to autonomously propose, evaluate, and optimize HDTwins. Specifically. LLMs iteratively generate novel model specifications, while offline tools are employed to optimize emitted parameters. Correspondingly, proposed models are evaluated and evolved based on targeted feedback, enabling the discovery of increasingly effective hybrid models. Our empirical results reveal that HDTwinGen produces generalizable, sample-efficient, and evolvable models, significantly advancing DT's' efficacy in real-world applications.
Introduction
Digital Twins (DTs) are computational models that accurately simulate the states and temporal dynamics of real-world systems [1, 2]. They particularly useful in modeling dynamical systems, which consist of multiple interdependent components that evolve over time [3, 4]. Take, for example, the epidemiological dynamics of contagious disease containing various components, including infection rates, recovery rates, population movement, and intervention strategies. DT's can integrate these factors to simulate future outcomes (e.g. predict disease spread), understand system changes (e.g. examining shifts in disease dynamics for varying demographics). and evaluate the impact of control measures (e.g. curb disease transmission) 6).
Desiderata. notable differentiator between DTss and general machine learning (ML.) models is the emphasis generalization. DTs are designed simulate completely unseen scenarios or interventions at inference time. Therefore, a crucial consideration [P1] out-of-distribution
'Equal contributions; authors listed in randomized order.
"odde available at https://gtthbb cwt/ http//mmoottHHDttiiniooo
38th Conference on Neural Information Processing Systems (NcurIPS 2024).

---

generalization: the ability to generalize to state-action distributions beyond those observed during training. This challenge is often compounded by the scareity of observational data available to accu- rately learn dynamics, highlighting the importance of > [P2] sample-efficient learning. Additionally, the model should be [P3] evolvable: capable of efficiently adapting (i.e. with minimal retraining) to changes in the underlying system dynamics. This is particularly crucial in healthcare domains, such as cpidemiological modeling and treatment planning, where DTss are regularly updated to reflect fundamental changes in discase transmission patterns (caused by viral mutations, vaccination coverage) or evolving drug resistance mechanisms, often with minimal additional data of emergent dynamics [6, 7].
Existing approaches for creating DTss primarily utilize two approaches: mechanistic models or ML- based neural models. Mechanistic models, denoted as fmech. are closed-form equations grounded in domain knowledge such as biological or physical principles. They offer high accuracy and generalization given sufficient domain understanding but are limited in their ability to model systems where scientific knowledge is incomplete [8, 9]. Of related note, techniques have been introduced to discover governing equations directly from data, but face challenges in scaling to more complex problem settings [10, 11]. Conversely, neural approaches, fncurat, leverage neural networks (NN) to learn DTs directly from data, often requiring minimal knowledge [12, 13, 14, 15]. Such models are effective given sufficient training data that provides adequate coverage of state- actiinn distributions, but struggle in data- data-scarce settings and are difficult to evolve to reflect changing conditions due to their overparameterized, monolithic nature.
Key considerations. Informed by this context, Hybrid Digial Twins (HDTwins) combine the strengths of both approaches through compositions of neural and mechanistic components, i.e. f = Smech fneural- Here, Smech symbolically incorporates domain- grounded priors, improving generalization and regularization while simplifying the complexity of patterns that have to be learned by the neural component. In other terms, fneural complements the mechanistic component by modeling complex temporal patterns in regions where the mechanistic model might be oversimplified or incomplete. Consequently, HDTwins can more accurately and robustly capture system dynamics, particularly in settings with (limited) empirical data and (partial) domain knowledge.
Conceptually, hybrid modeling involves two stages: model specification, determining the model structure (e.g. neural architecture, symbolic equations), and model parameterization, estimating model parameters (e.g. neural weights, coefficients). This process, with model specification in particular, has traditionally relied heavily on human expertise to craft [16, 17, 18, 19]. In this work, we investigate the feasibility of automatically designing hybrid models with minimal expert involvement, which would significantly enhance the efficiency and scalability of model development. This task is challenging, as it requires searching for optimal specification and corresponding parameters within a vast combinatorial model space [20, 21J. To address this, we introduce HDTwinGen, novel evolutionary framework that autonomously and efficiently designs HDTwins. At a high level, our method represents hybrid model specifications in code and leverages large language models (LLMs) for their domain knowledge, contextual understanding, and learning capabilities to propose symbolically represented models and search the model space [22, 23, 24). This is coupled with offline optimization tools to empirically estimate model parameters from training data. More specifically, HDTwinGen utilizes two LLM agents: the modeling agent, whose task is to generate novel model specifications. and the evaluation agent, which analyzes performance and provides targeted recommendations for improvement. Through multiple iterations, HDTWwinGen efficiently evolves better performing hybrid models with informed modifications.
models
Contributions: ! Conceptually, we present the first work in automated hybrid model design, jointly optimizing model specification and parameterization of hybrid digital twins. 2 Technically, introduce HDTwinGenn a novel evolutionary framework employing LLMs and offline optimization tools to propose, evaluate, and iteratively enhance hybrid models. 3 Empirically, we demonstrate that our method learns more accurate D'Ts, achieving better out-of-distributtion generalization, sample-efficient learning, and increased flexibility for modular evolvability.
2 Digital Twins of Dynamical Systems
A dynamical system S = (X, Φ) is tuple of its dx-dimensional state space X C Rxx, an (optional) dimiinnnsol action space U C Rdd, and dynamics model Φ. The state at time

---

E T C R4 is represented as a vector, x(t) E X' and similarly the action taken is represented as vector u(t) U. The continuous-time dynamics of the system can be described by dx(t) (xx(t), u(t), I), where + X x U x T + X. We optionally consider the existence of some policy T: :: X + P(U) that acts on the system by mapping a state sr()) to a distribution over actions u(t).
Digital Twins. Digital twins (DTs) aim to approximate & X x u x T + X' using a computational model (o) E F learned from data. Here, we use 0 E θ to denote the specification of the model (e.g. linear) and (()) (2()) to indicate the set of parameters specified by 8. Additionally. F, O, and (2) are the spaces of all possible models, specifications, and parameters, respectively. Next, we outline the key desiderata for DT:
[PI] Generalization to unseen state-action distributions. As DTs are required to simulate vary- ing conditions, they should extrapolate to aate-action distributions not observed during training time. Formally, the generalization error ((z(t), ((()))) ((t))))))))) ( [C(So (x(t), u(t)), u(t)) should be minimized, where C is some loss function, and PooD represents the outtof-distribution scenario. [P2] Sample-efficient learning. Given the often limited availability of real-world data, DT's should learn robustly from minimal empirical data. In other words, they must have good sample complexity, achieving the desired level of generalization with a limited number of observations [25].
[P3] Evolvability. Dynamical systems are, by nature, non-stationary and evolve over time [26, 27]. From a modeling perspective, the DT should be easily evolved to reflect changing underlying dynamics, minimizing the need for additional data or expensive model re-development, i.e. 0 and (0) should be easily adjustable to reflect changing system dynamics.
For the purpose of model learning, we assume access to offline dataset containing NN E N+ trajec- tories, where the measurements the systems are made at discrete time points [7] = [2,,,, f2. ..... T7. This dataset, D = {{(=(() (t), u (t), ((m) (t)) | t E [T])]== contains state-action trajectories sampled regularly over time, where (t) = x (m) (t + At) represents the subsequent state.
3 Hybrid Digital Twins
HDTwin. A Hybrid Digital Twin is composition of mechanistic and neural components, represented as fo, = fmech o facural [18, 28). This class of hybrid models offers several advantages that align with our desiderata. The mechanistic component allows partial knowledge to be encoded through its symbolic form, which, while not sufficient alone to accurately predict underlying dynamics, is complemented by the neural components that learn from available data. This combination aids in generalization ([P1)), especially moving beyond conditions observed in training, and improves sample complexity ((P2). Furthermore, the mechanistic component can be quickly and easily updated with new parameters due to its simpler, lower-dimensional structure, allowing the overall model adapt efficiently to remain accurate in changing conditions ((P3). In this work, we focus on additive compositions, fa, = f + fueurat. as they are more interpretable. Additionally, it enables individual contributions of mechanistic and neural components to easily disentangled and simplifies the optimization to allow gradient-based methods [16]. Nonetheless, we encourage future works to investigate alternative composition strategies (e.g. branching composition) to develop more advanced HDTwwis [29].
Learning the hybrid model can be decomposed into steps: (1) model specification, learning the structure, 0 E θ, of the dynamics function that describes how the system evolves over time; and (2) model parameterization, which estimates the specific values of parameters w0) 22))) for a given specification 0. For instance, the logistic-growit model specifies a structure for population growth, while parameterization involves estimating the growth rate and carrying capacity." 3 More generally, this learning problem can be mathematically formulated as a bilevel optimization problem:
8* arg min Comer w" (0)), where dEE
w- arg min (0, (+)) wE
(1)
Here, the upper-level problem involves finding the optimal specification 0* that minimizes the outer objective Louter, while the lower-level problem involves finding the optimal parameters w'(() for
given specification 0 that minimizes the inner objective function Linne- To be more concrete, the
TThis model relates population size N. growth rate r, and carrying capacity K, aN/oc rNN N/xx) [30].
3

---

1) Modelling generales Seucdy 1.N 16 pia-11, (ia-- gonnttt sporifteation
Modeling
User
grontest
Contect
2) Optimioon mocial paramutetiaation
H1a-1)
Sountert
Modeling LLNmadel
206 argeilfy Aluwo
Fraluation Pardlackk
Mudilliu Contect
Set of Topak HDTwin
HDTwin Jaw
Evaination agent geverates 1110) LLM (chontter plat)
NI. Foodback
Update population
Evaluate HDTwia
Evaliaation
Soontest
gist Wo tt s
Elfewar Daa)
Compute amrnage liwe
LLMeval
Agent
Un ((0 HDTwim
plo
-pp--ll
6 diaaiey Pan))
Compure compooanti
Conrect
U
Figure 1: HDTwinGen: evolutionary framework. The process begins with user-provided modeling context gromiext and D = {Duain }. 1) In iteration 9. the modeling agent generates model specification Python program ,0 2) Parameters are optimized using the offline optimization tool to yield fo.-~ (0). 3) The HDTwin is evaluated based on model loss v and component-wise loss 8. Subsequently, the model pool pla) is updated with top-K models. 4) The evaluation agent provides targeted feedback for model improvement H by analyzing models in P(s) using performance metries requirements outlined in Scontext This iterative loop repeats for G iterations.
outer objective measures the generalization performance, empirically measured on the validation set Lal while the inner objective measures the fitting error, as evaluated the training Liain-
Combinatorial search space. The space of possible specifications θ (e.g. different networks, functional forms) is discrete and combinatorially large, while (2()) represents the continuous space of parameters to be optimized. Selecting the optimal 0, (() thus involves searching through a vast combinatorial space. Performing this search through traditional means, such as genetic programming [21] or evolutionary algorithms [20], is computationally challenging, time-consuming, and often technically infeasible. To the best of our knowledge, automatic HDTwin development, where we incorporate LLMs (combined with offline optimization tools) to automatically optimize both the specification and the parameterization of hybrid models.
work is the first to address the problem of
4 HDTwinGen: Automatic Design of HDTwins
Human experts craft models by making strategic design decisions based on their domain knowledge, starting with sensible initial model specification and performing intelligent modifications based on empirical evaluations. Our key insight is that LLMs can effectively emulate these capabilities to efficiently navigate the search space in Equation (1) and autonomously design HDTwins. More specifically, our method utilizes LLMs for three major purposes: source of domain knowledge, where LLMs inject consistcnnnsstent knowledge into the model specification, particularly through the symbolic representation fmech efficient search, by making intelligent modifications to the specification to converge more efficiently on the optimal hypothesis; and contextual understand- ing, enabling the algorithm to incorporate task-specific context and targeted feedback for model improvement [22, 23, 24].
We operationalize this insight through HDTwinGen, evolutionary algorithm that iteratively evolves population of candidate solutions to automatically search for the best HDTwin. Our approach employs framework comprising three key elements: (I) human experts provide an initial system description, modeling objectives, and requirements as a structured prompt; (2) modeling agent proposes new model specifications, optimizes their parameters dataset, and collects validation performance metrics; (3) an evaluation agent assesses the proposed models using both data-driven performance metries and qualitative evaluations against expert-defined objectives and requirements. The agents communicate using natural language and a custom code format representing the HDTwin model, facilitating autonomous and iterative model enhancement. An overview of our method is presented in Figure 1, with pseudocode in Appendix E.
Overview.
training
Initial prompt design. The optimization process begins with human expert providing a structured prompt, referred to as the modeling context S gooneet This modeling context outlines the system description, modeling objectives L, and requirements R:

---

1. The system description semantically describes the system, including state and action variables, giving the algorithm the contextual understanding necessary for informed model development. The modeling objective specifies quantitative performance requirements via a metric L. The modeling requirements R are qualitative and described in natural language, detailing aspects such as interpretability (e.g. fully mechanistic or hybrid model) and additional scientific knowledge (e.g. a log-linear relationship between variables).
In practice, R. can incorporate various requirements, allowing for the design of both purely mechanis- tic and hybrid models, a flexibility that we demonstrate experimentally. The model is represented in Python, where purely mechanistic specifications are represented in native Python and neural components are represented using PyTorch [31]. Moreover, gemmiexi includes a skeleton code to guide the synthesis of executable code in predetermined format. For illustrative purposes, an example of goontext is provided in Appendix E.4.
Evolutionary optimization overview. Given Scontext as input, HDTwinGen performs G' iterations of optimization, where G E N+ The population of proposed HDTwins at iteration g is represented as p(s), Each iteration creates new candidate model based on previously created models in p(s) and feedback. Only the top K models are retained after each iteration, except when g < K, in which case all generated models are kept, i.e. max [p(0)] K. Each model in P(a) is characterized by a tuple containing its model specification (represented symbolically through code) and validation metrics. After completing G iterations, the model with the best validation performance in p(0) is selected as the final model.
4.1 Modeling Agent
The goal of the modeling step is to propose novel HDTwins based on previously proposed models and feedback from the evaluation agent. Specifically, on the g-th iteration, the modeling agent takes as input * pls-1); the set of top-K previously generated models; H I); the most recent feedback produced by the evaluation agent (where the initial step, g = both are empty, i.e., (0) Ø. p(0) Ø); and + geomlext, the modeling context. The modeling agent generates a model specification 0 using a predefined code format (i.e. skeleton code). By observing multiple previously best-performing models and their performances, the modeling agent can exploit this context as a rich form in context-learning and evolve improved specifications in subsequent generations [22]. Each generated specification emits its corresponding parameters, w()) are fitted to the training set Duain. More formally, we represent this generative procedure as fo. ~ LLMmodel (H, p(9), Scomlext).
Proposing HDTwins.
Model specification. To generate model specifications, the modeling agent decomposes the system into a set of components, with each component describing the dynamics of specific state variable. In other words, for system with dx state variables, there will be dx components. Each component is characterized by its own set of inputs and a unique dynamics function that describes the dynamies of its associated state variable over time. This modular representation enables independent analysis and optimization of individual components. cases where R specifies purely mechanistic equations, the component dynamics are entirely defined using closed-form equations. Conversely, in hybrid model, the mechanistic equation can be augmented with a neural network (implemented in PyTorch to model residuals (i.e. in additive fashion). The choice between mechanistic and hybrid models is left to the user, balancing the trade-off between transparency and predictive performance. Concretely, the specification step involves 'filling in' the skeleton code with a detailed body of code, specifying the decomposition, and delineating each component's 's dynamics function as a separate code structure (for generated HDTwin example, see Appendix l).
Model optimization. The generated specification emits (0), which are treated as placeholder values, and are then optimized against the training dataset. Specifically, we optimize the mean squared error for the parameters that minimize this loss as w (0) = arg min [(fo Duain). In this work, we consider (8) to be continuous variables, and as such, we optimize 0 by stochastic gradient descent, using the Adam optimizer [32]. However, we note other optimization algorithms, such as blxck-box optimizers, could also be used (for more details, see Appendix F, Equation (5)). The parameter optimization step then yields the complete model, (0)
Quantitative evaluation. For each generated model, we evaluate them quantitatively. Specif- ically, we collect the validation mean squared error loss per component, which we denote as
5

---

....,,..... ,,, (Appendix F, Equation (6)). We also compute the validation loss of the 5 3- =- [01,02,- overall model as well as v = (() Daa)). Finally, the generated model and its validation losses are included in a tuple and added to the top- K models ple) + p(s-1) @ (fo,w-(0), 6, v), where pla) automatically removes the lowest performing models, and also only adds new model to pla) if is unique. We highlight that we consider the top-K models only to apply selection pressure, such that only the best-performing models are considered when generating the next HDTwin [33].
4.1.1 Evaluation Agent
The goal of the evaluation step is to reflect on the current set of top- K models, P(s) against requirements R and provide actionable and detailed feedback to the modeling agent for model improvement: H(9) ~ LLMesai P(s)). We note that H'9 is provided in natural language and can be viewed as dense feedback signal, a notable distinction from traditional learning methods, where feedback often takes the form of simple scalar values, such as loss gradients or rewards. Leveraging natural language feedback allows the agent to (I) engage in comparative analysis, identifying effective specifications in po()) contributing higher performance and discerning patterns common in less effective models, informing its suggestions for further model improvement; (2) qualitatively evaluate models against qualitative requirements R -eveerraggn the LLM's capacity to reason about proposed HDTwins to reflect these requirements via model improvement feedback.
Model evaluation.
Enhancing search. By providing rich feedback improve model specification, the evaluation and modeling agent collaborate to efficiently evolve high-performing models. Empirically, in Appendix J, we observe that the evaluation agent provides targeted and specific feedback, including component- specific suggestions, proposing alternative decompositions, removing parameters, or introducing non-linear terms. It is noteworthy that the feedback H()), expressed flexibly in natural language, could easily be further enriched through direct human feedback. We demonstrate this human-in- thoploop capability by including expert feedback during the optimization process through (9)) and observed that it was integrated into newly generated HDTwins. Though further investigation is beyond the scope of this work, this demonstration highlights promising avenues for augmenting human-machine collaboration in the autonomous design of DTTs.
5 Related Works
For an extended related work, refer to Appendix B. Our work focuses on autonomously learning DTs from data, with several relevant research strands:
Neural sequence models. ML approaches commonly address learning system dynamics se- quential modeling problem. In these settings, fo. are typically black-box models, where 0 E θ the NN architecture and w() are its weights. Early models like Hidden Markov Models [34] and Kalman filters [35] made simplifying Markovian and lincarity assumptions, later extended to nonlinear settings [36, 37J. Subsequent models, including recurrent neural networks [38]. along with their advanced variants [39, 40, 41], introduced the capability to model longer-term depen- dencies. More recent advancements include attention mechanisms [42] and Transformer models [43], significantly improving the handling of long-term dependencies in sequence data. Another line of work, Neural Ordinary Differential Equations (NODE) [14, 44, 45], interprets neural network operations differential equations. These methods have found utility in modeling range of complex systems [46, 47, 48, 49). While deep sequence models are proficient at capturing complex dynamics, they are heavily reliant on training data for generalization (PP1, P2]), and their monolithic and overparameterized structures limit evolvability ([P3)).
Mechanistic (discovery) models. Beyond purely neural approaches, another line of work aims to discover a system's s governing equations directly from data. Here 0 are closed-form equations and (() are their parameters. These include symbolic regression techniques [10]. Eureqa [50], SINDY [11], D-CODE [51, 52], among others [52, 53] that search for @ and w(f) from data. These techniques struggle to scale digher-dimensiooal settings and rely experts perform variable selection and define the function set and primitives available to the search algorithms.
Hybrid models. Recent efforts have also created hybrid models by integrating physical laws with neural models. Physics-informed neural networks [15, 54], and methods including Hamiltonian Neural Networks [55], Lagrangian Neural Networks [56] integrate structural priors of physical

---

systems to improve generalization. These techniques introduce specialized mechanisms to incorporate precisely known physical principles. Additionally, [57] integrates prior ODE/PDE knowledge into hybrid model, using specialized regularization to penalize the neural component's information content. [58, 59] consider settings where an expert equation is known, but equation variables are latent and unobserved. Correspondinglyy they employ two sets of latent variables: one governed by expert equations and another linked to neural components. [60] performs data augmentation by sampling out-of-distribution trajectories from expert models. While existing approaches rely on expert models to perform the hybrid model design, HDTwinGen an automated approach to jointly optimize hybrid model specification and its parameters.
6 Experiments and Evaluation
In this section, we evaluate HDTwinGen and verify that it significantly outperforms state-of-ther afrt methods in modeling system dynamics over time from observed dataset and corresponding system description." 4
Benchmark datasets. We evaluate against six real-world complex system datasets; where each dataset is either a real-world dataset or has been sampled from an accurate simulator designed by hu- man experts. Three are derived from state-of-the-ar biomedical Pharmacokinetic-PP aaapeookineii--Paaaaaadyyyeee (PKPD) model of lung cancer tumor growth, used to simulate the combined effects of chemotherapy and radiotherapy in lung cancer [61] (Equation -2))-this has been extensively used by other works [62, 63, 64]. Here we use this bio-mathematical lung cancer model to create three variations of lung cancer under the effect of no treatments (Lung Cancer), chemotherapy only (Lung Cancer (with Chemo.)), and chemotherapy combined with radiotherapy (Lung Cancer (with Chemo. & Radio.)). We also compare against an accurate and complex COVID-19 epidemic ggent-based simulator (COVID-19) [65], which is capable of modeling ao--harmmaeuutcall interventions, such as physical distancing during lockdown. Furthermore, we compare against ecological model of a microcosm of algae, flagellate, and rotifer populations (Plankton
Pharmacodynamic
(erpecssmmerrlliattinn an experimental three-species prey-predator system [66]. Moreover, we also compare against real-world dataset of hare and lynx populations (Hare-Lynx), replicating predator-prey dynamics [67]. We detail all benchmark datasets details in Appendix C.
Evaluation Metrics. We employ mean squared error (MSE) to evaluate the benchmark methods on a held-out test dataset of state-action trajectories, denoted Drest: using the loss defined in Equation (5) and report this as TMSE. Each metric is averaged over ten runs with different random seeds, and we present these averages along with their 95% confidence intervals, further detailed in Appendix G.
Benchmark methods. To assess whether HDTwinGen is state-of-the-an we compare it with the most competitive and popular neural network models, which, when modeling the dynamics of system over time, becomes form of ODE model, that is a neural ODE [14] with action inputs (DyNODE) [68]. Moreover, we also compare against recurrent neural network (RNN) [69] and state-of-the-art transformer (Transformer) [64]. We also compare against mechanistic dynamical equations derived from equation discovery methods for ODEs, including Genetic Programming (GP) [10] and Sparse Identification of Nonlinear Dynamics (SINDy) [11]. Lastly, we compare against hybrid model (APHYNITY) that integrates prior knowledge in the form of ODEs into hybrid models, while penalizing the information content from the neural component [57]. Moreover, we compare against the ablations of our method, of the zero-shot generated HDTwin (ZeroShot) and this model with subsequently optimized parameters (ZeroOptim). We provide method implementation, hyperparameter, and experimental details in Appendix D.
7
Main Results
We evaluated all our benchmark methods across all our datasets tabulated in Table 1. HDTTwinGen models the system the most accurately, achieving the lowest test prediction mean squared error on the held-out test dataset of state-action trajectories. In the interest of space, we include additional experimental evaluations in the appendix. Specifically, we also evaluate HDTwinGen performance on a suite of synthetically and procedurally generated benchmarks (Appendix H.9); > comparisons
'Our implementation is available at /tttps///gtthbb.oo. http:///gg.hat...oooooooottttttt We also provide a wider
lab code repository https: //github com hcn///aandrrrehaaaaaa///iiaaaa
7

---

Table 1: Benchmark method performance. Reporting the test prediction MSE (TMSE) of the produced system models on held-out test datasets across all benchmark datasets. HDTwinGen achieves the lowest test prediction error. The results are averaged over ten random seeds, with + indicating 95% confidence intervals.
Lung Cancer
Ling Canser (with Chenio))
Lung Cancer (with Chetio & Kidic.)
Hart Lynx
Plankion Microcosnn
CONTES
Medeed
Tuas
Tarse
Turse
Tarza
Tiess
Tarae
SINEy
79
11.
IXVEILS13
355:x227-11
BIOOESEU
GP
171-899
S14+
DSNODE
15,0E533
139210
742226
LITOM0R:3 3.04440
13545.66
3.71 30u-03
0.038 +0.400
Trandommer
7486:100
0:2482000018
716a4425
APHYNUIY
M00631.37
$1:0:81.3
1.21e4036:
331::12.66
5.81e403944
ZenoOplin
2162:72
31,21:45
6.05.37.9
353.10
0,0132::000166
7.88:000011
HDTTwinGen
0,6889 1UN453
0,131
201 38.3
2,51e-80 1.22e-66
1.22
Table 2: Out of distribution shifts. On varia-
tion of the Lung Cancer (with Chemo. & HDTwinGen is more robust OOD shifts in un- seen state-action distributions,
Radio.),
100
10
10F
10 10F
10
Nardon of Incedor..
Lung Cancer
(with Chemo. & Radio.)
Figure 2: Sample efficiency. Analyzing perfor- mance function of the number of training trajectories in the Lung Cancer (with Chemo. & Radio.) dataset. We observe that HDTwinGen achieves the lowest test prediction error, even in the very challenging low data regime. This high- lights the role of priors embedded in HDTwin in sample-efficient generalization.
Method
IID Tuss 4
0OD Tusss
DyNODE
0.011320.0121
1.75=0.769
SINDY
0.302:10226
5.922.55
RNN
448e=04442.112ewl
Transformer
0.0262+00.0514
199==44:-2.8e440)
ZeroShot
4.55=03++1.4544000 M.91
ZeroOptim
3.49:1000364
4.844:5.17
HDTwinGen
0.00872:10.0187
0.0846:10.0891
against domain-specing baselines (Appendix H.8) and various ablation experiments, including ablation of LLM hyperparameters, prompt design, and algorithm settings (Appendices H.5 to H.7).
7.1 Insight Experiments
This section provides in-depth analysis of HDTwinGen''s effectiveness related to its benchmark counterparts. Specifically, we examine the core desiderata for an effective DT described in Section 2:
[P1] -utto--distribution generalization, [P2] sample-efficient learning, and [P3] evolvability.
[P1] Can an HDTwwin generalize to adapt the Lung Cancer (with Chemo. & Radio.) produce training dataset of states in range that is outside those observed in the test set over all trajectories (Appendix H. We tabulate this in Table 2. Empirically, we find that HDTwinGen is more robust to out-of auttof-distribution shifts than existing methods, benefiting from explicit decomposition and robust hybrid models. Notably, the neural network method DyNODE shows the largest relative error increase from IID to OOD by two orders of magnitude, while the mechanistic method SINDy exhibits a smaller increase by only one order of magnitude. This demonstrates the importance of hybrid models that leverage both neural and mechanistic components to enhance generalization performance under distribution shifts.
shifts? To explore
shifts,
[P2] Can HDTwinGen improve sample-efficiency in model learning? To explore the low data settings, we re-ran all benchmark methods with fewer samples in their training dataset on the Lung Cancer (with Chemo. & Radio.) dataset. We plot this in Figure 2. Empirically, observe that HDTwinGen can achieve lower performance errors, especially in lower-sample regimes.
[P3] Can HDTwinGen evolve its modular HDTwin to fit the system? We analyze this from an empirical point of view to determine if HDTwinGen can correctly evolve the generated HDTwin and reduce its prediction eтror over subsequent generations. We observe that HDTwinGen can indeed understand, reason, and iteratively evolve the generated code representation of the HDTwin to incorporate better fitting HDTwin, as observed in Figure 3. In particular, the annotated results demonstrate that HDTwinGen effectively refines the hybrid model by strategically adjusting its neural and mechanistic components (in a fashion akin to human experts), leading to significant improvements
8

---

cmbirnd Phengy diie herepaenntl Jony draamoderryyy cnnettzzio
Tope HDTxt
HDTTwi Geeettel Sap
Validation MSE
Inital whii bet
100
Whi biack mool lepiaic
black
logistic
Merkon
Merun orumatio chamosherrpy effecr, eggendd ferzion radurturayy effet, resitanee deckkgrrr, anl arsplite shok bes MLL roadaal...
effes, egrond forctin for radiethenpy etfec, reoeazee developmer, Mack bos MLP reskkalss.
White
10-1
dynamtc carying canacity and donoy tin, non-itezar insoccioon omm, spebeed back hove ML.P Sernesatte
4
6
8
10
12
14
Number of Generations
Figure 3: HDTwinGen effectively evolves HDTwin. Validation MSE the HDTwin generated in each iteration, showing the Pareto-front of the best generated HDTwin (Top-1 HDTwin), and the generated HDTwin per generation step-adee addiiioiio with a few of the HDTwins labeled with their model descriptions. HDTwinGen can efficiently understand, modify, and hence evolve the HDTwin to achieve a better-fitting model (Appendix H.4).
Exposed (# Individuals)
Sudden Intervention of Lockdown
100000
50000
Observed and trained on
Unobserved to predict
0
5
10
15
20 (Days)
25
30
35
40
True
HDTTwin (with intervention)
SINDy
DyNODE
Figure 4: COVID-19 unobserved intervention. The symbolic code-based representation of HDTwin can be easily adapted to unobserved interventions through targeted adjustments of parameters.
accuracy and robustness. This iterative evolution process demonstrates HDTWinGGen's is ability to adapt and optimize its modular components.
Can HDTwinGen Understand and Modify Its HDTwin? We investigate whether large language model (LLM) agents can take an optimized high-dimensional twin (HDTwin) from an existing benchmark dataset and adapt model an unobserved intervention that is not present in the training data. We note that this intervention emulates scenarios where the dynamics of the underlying system changes. We affirmatively answer this question by constructing a scenario where our COVID-19 simulator incorporates unobserved intervention of lockdown policy, which reduces physical interactions between individuals (Appendix H.2). As demonstrated in Figure 4, we observe that the code-model representation of the HDTwin can be (1) understood by the modeling agent LIM and (2) adapted its parameters to accurately model and reflect this intervention. We find that HDTwinGen is the only method capable of changing the overall functional behavior by modifying single parameter in the model; in contrast, all other existing data-driven methods require a dataset of state-action trajectories under the new dynamics introduced by this intervention.
Ablation Studies. We conducted ablation studies on HDTwinGen and found several key insights. First, retaining the top-K models within the LLM context leads to improved model generation (Ap- pendix H.5). Additionally, HDTwinGen compatible with various LLMs and different temperature settings (Appendix H.6). It also benefits from including textual descriptions of the variables to be modeled as prior information (Appendix H.7). Finally, HDTwinGen can be specifically instructed to generate mechanistic white-box models if desired (Appendix 10).

---

8 Limitations and Discussions
In summary, this work addresses the problem of learning digital twins for continuous-time dynamical systems. After establishing clear learning objectives and key requirements, we introduce Hybrid Digital Twins (HDTwins)) promising approach that combines mechanistie understanding with neural architectures. HDTwins encode domain knowledge symbolically while leveraging neural networks for enhanced expressiveness. Conventional hybrid models, however, rely heavily on expert specification with learning limited to parameter optimization, constraining their scalability and applicability. To overcome these limitations, we propose a novel approach to automatically specify and parameterize HDTwins through HDTwinGen, an evolutionary framework that leverages LLMs to iteratively search for and optimize high-performing hybrid twins. Our empirical results demonstrate that evolved HD)Twins consistently outperform existing approaches across multiple criteria, exhibiting superior dii-stiiittittno generalization, enhanced sample efficiency, and improved modular evolvability.
Limitations. While our results are promising, several important limitations remain. HD'winGee'' efficacy depends critically on human experts providing initial system specifications and on the underlying LLM's 's domain knowledge and model generation capabilities. Our current implementation focuses exclusively on continuous-ting titooouss-tteeeeee systems, which, although broadly applicable, represent only subset of real-world systems. Future work could extend our approach through human-in-thee the-loop feedback mechanisms, integration with external tools, and expansion to broader system classes.
Ethical implications. We acknowledge the risk of bias transmission from the black-box LLMs into the evolved models. While our hybrid approach enables greater expert scrutiny through its human-interpretable components, we strongly recommend comprehensive evaluation of evolved models for fairness, bias, and privacy concerns before deployment in sensitive applications.
Acknowledgments and Disclosure of Funding
We thank the anonymous reviewers, and program chairs, members of the van der Schaar lab, and Andrew Rashbass for many insightful comments and suggestions. TL and SH would like to acknowledge and thank AstraZeneca for their sponsorship and support. This work was supported by Microsoft's Accelerate Foundation Models Academic Research initiative.
References
[1] Tao, J. Cheng, Q. Qi, M. Zhang, H. Zhang, and Sui, 'Digital twin-driven product design, manufacturing and service with big data," The International Journal of Advanced Manufacturing Technology, vol. 94, pp. 3563-3576, 2018
[2] Corral-Acero Margara, M. Marciniak, C. Rodero, Lonearic, Y. Feng. A. Gilbert, F. Fernandes, H. A. Bukhari, A. Wajdan et al., "The 'digital twin't to enable the vision of precision cardiology." European heart journal, vol. 41, no. 48, pp. 4556-4564, 2020.
[3] H. A. Simon, The sciences of the artificial. MIT press. 1996.
[4] J. Ladyman, J. Lambert, and K. Wiesner, "What is complex system?" European Journal for Philosophy of Science, vol. 3, pp. 33-67, 2013.
[5] Q. Qi and Tao, "Digital twin and big data towards smart manufacturing and industry 4.0: 360 degree comparison." leee Access, vol. 6, pp. 3585-3593, 2018.
[6] V. Iranzo and S. Pérez-Gonzálezz. "Epidemioolgical models and covid-19; a comparative view," History and Philosophy of the Life Sciences, vol. 43, no. 3, p. 104, 2021.
[7] I. Bozic, J. G. Reiter, B. Allen, T. Antal, K. Chatterjee, P. Shah, Y. S. Moon, A. Yaqubie, N. Kelly, D. T. Le et al., "Evolutionary dynamics of cancer in response to targeted combination therapy." elife, vol. 2, p. e00747, 2013.
[8] R. Rosen, G. Von Wichert, G. Lo, and K. D. Bettenhausen, "About the importance of autonomy and digital twins for the future of manufacturing." sfac-papersonline, vol. 48, no. 3, pp. 567-572, 2015.
10

---

[9] T. Erol, A. F. Mendi, and D. Doğan, "The digital twin revolution in healthcare," in 2020 4th international symposium on multidisciplinary studies and innovative technologies (ISMSIT). IEFE, 2020, pp. 1-7.
[10] J. R. Koza, "Genetic programming as a means for programming computers by natural selection," Statistics and computing. vol. 4, pp. 87-112, 1994.
[11] S. L. Brunton, J. Proctor, and N. Kutz, "Discovering governing equations from data by sparse identification of nonlinear dynamical systems," Proceedings of the national academy of sciences, vol. 113, no. 15, pp. 3932-3937, 2016.
[12] D. Ha and Schmidhuber, "Recurrent world models facilitate policy evolution," Advances in neural information processing systems, vol. 31, 2018.
[13] Yoon, D. Jarrett, and M. Van der Schaar, "Time-series generative adversarial networks," Advances in neural information processing systems, vol. 32, 2019.
[14] R. T. Chen, Y. Rubanova, J. Bettencourt, and D. K. Duvenaud, "Neural ordinary differential equations," Advances in neural information processing systems, vol. 31, 2018.
[15] M. Raissi, P. Perdikaris, and G. E. Karniadakis, "Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations." " Journal of Computational physics, vol. 378, pp. 686-707, 2019.
[16] L. Faure, B. Mollet, W. Liebermeister, and J.-L. Faulon, "A neural-mechanistic hybrid approach improving the predictive power of genome-scale metabolic models," Nature Communicutions, vol. 14, no. 1, p. 4669, 2023.
[17] Pinto, R. Ramos, R. S. Costa, and R. Oliveira, "A general hybrid modeling framework for systems biology applications: Combining mechanistic knowledge with deep neural networks under the sbml standard," AI, vol. 4, no. 1, pp. 303-318, 2023.
[18] Wang, Z. Zhu, W. Liang, L. Liao, and J. Wan. "Hybrid mechanistic and neural network modeling of nuclear reactors," Energy, vol. 282, p. 128931, 2023.
[19] R. Cheng, A. Verma, G. Orosz, S. Chaudhuri. Y. Yue, and J. Burdick, "Control regularization for reduced variance reinforcement learning," " in International Conference on Machine Learning. PMLR, 2019, pp. 1141-1150.
[20] E. Real, S. Moore, A. Selle, S. Saxena, Y. L. Suematsu, J. Tan, Q. V. Le, and A. Kurakin, "Large-scale evolution of image classifiers," in International conference on machine learning. PMLR, 2017, pp. 2902-2911.
[21] T. N. Mundhenk, M. Landajuela, R. Glatt, C. P. Santiago, D. M. Faissol, and B. K. Petersen, "Symbolic regression via neural-guided genetic programming population seeding," in Proceed- ings of the 35th International Conference on Neural Information Processing Systems, 2021, pp. 24 24912-24 923.
[22] T. Brown, B. Mann, N. Ryder, M. Subbiah, J. D. Kaplan, P. Dhariwal, A. Neelakantan, P. Shyam, G. Sastry, A. Askell el al., "Language models are few-shot learners," Advances in neural information processing systems, vol. 33, pp. 1877-1901, 2020.
[23] J. Wei, X. Wang, D. Schuurmans, M. Bosma, F. Xia, E. Chi, Q. Le, D. Zhou et al., "Chain-of- thought prompting elicits reasoning in large language models," Advances in neural information processing systems, vol. 35, pp. 24 824 24 24 837, 2022.
[24] A. Chowdhery, S. Narang, J. Devlin, M. Bosma, G. Mishra, A. Roberts, P. Barham, H. W. Chung, C. Sutton, S. Gehrmann et al., "Palm: Scaling language modeling with pathways," Journal of Machine Learning Research, vol. 24, no. 240, pp. 1-113, 2023.
[25] M. J. Kearns and U. Vazirani, An introduction to computational learning theory. MIT press, 1994.
[26] H. A. Simon, "The architecture of complexity," Proceedings of the American philosophical society, vol. 106, no. 6, pp. 467-482, 1962.

---

[27] T. L. Rogers, B. J. Johnson, and S. B. Munch, "Chaos is not rare in natural ecosystems," Nature Ecology & Evolution, vol. 6, no. 8, pp. 1105-1111, 2022.
(28) M. Sokolov, M. von Stosch, H. Narayanan, F. Feidl, and A. Butté, "Hybrid modeling key enabler towards realizing digital twins in biopharma?" Current Opinion in Chemical Engineering, vol. 34, 100715, 15, 2021.
[29] S. Chaudhuri, K. Ellis, O. Polozov, R. Singh, A. Solar-Lezama, Y. Yue et al.. "Neurosymbolic programming," Foundations and Trends@ in Programming Languages, vol. 7, no. pp. 158- 243, 2021.
[30] A. Tsoularis and J. Wallace, "Analysis of logistic growth models," Mathematical biosciences, vol. 179, no. 1, pp. 21-55, 2002.
[31] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga et al., "Pytorch: An imperative style, high-performance deep learning library," Advances in neural information processing systems, vol. 32, 2019.
[32] D. P. Kingma and J. Ba, "Adam: A method for stochastic optimization," arXiv preprint 14Xiv:1412.6980 2,6900 2014.
[33] J. H. Holland, "Genetic algorithms," Scientific american, vol. 267, no. 1, pp. 66-73, 1992.
[34] L. R. Rabiner, "A tutorial on hidden markov models and selected applications in speech recognition," Proceedings of the IEEE, vol. 77, no. 2, pp. 257-286, 1989.
[35] R. Kalman, "A new approach to linear filtering and prediction problems," Trans. ASME, D. vol. 82, pp. 35 34, 1960.
[36] L. Li, Y. Zhao, D. Jiang, Y. Zhang, F. Wang, I. Gonzalez, E. Valentin, and H. Sahli, "Hybrid deep neural network-hidden in 2013 Humaine association conference on affective computing and intelligent interaction. IEEE, 2013, pp. 312-317.
hidden markov model (dnn-hmm) based speech emotion recognition,"
[37] R. G. Krishnan, U. Shalit, and D. Sontag, "Deep kalman filters," arXiv preprint 1rXiv:1511.05122,, 2015.
[38] J. L. Elman, "Finding structure in time." Cognitive science, vol. 14, no. 2, pp. 179-211, 1990.
[39] S. Hochreiter and J. Schmidhuber, "Long short-term memory." Neural computation, vol. no. 8, pp. 1735-1780, 1997.
[40] K. Cho, B. Van Merriënboer, C. Gulcehre, D. Bahdanau, F. Bougares, H. Schwenk, and Y. Bengio. "Learning phrase representations using rnn encoder-decoder for statistical machine translation," " arXiv preprint rrXiv:1406.1077 2014.
[41] I. Sutskever, O. Vinyals, and Q V. Le, "Sequence sequence learning with neural networks." Advances neural information processing systems, vol. 27, 2014.
[42] D. Bahdanau, K. Cho, and Y. Bengio, "Neural machine translation by jointly learning to align and translate," arXiv preprint arXiv: 1409:0430 2014.
[43] A. Vaswani, N. Shazeer, N. Parmar, Uszkoreit, L. Jones, A. N. Gomez, Ł. Kaiser, and I. Polosukhin, "Attention is all you need," Advances in neural information processing systems, vol. 30, 2017.
[44] E. Dupont, A. Doucet, and Y. W. Teh, "Augmented neural odes," Advances in neural information processing systems, vol. 32, 2019.
[45] S. I. Holt, Z. Qian, and M. van der Schaar, "Neural laplace: Learning diverse classes of differential equations in the laplace domain," in International Conference on Machine Learning. PMLR, 2022, pp. 8-832.
12

---

[46] M. A. Zaytar and C. El Amrani, "Sequence to sequence weather forecasting with long short-term m recurrent neural networks," International Journal of Computer Applications, vol. 143, no. 11, pp. 7-11, 2016.
[47] J. Devlin, M.-W. Chang, K. Lee, and K. Toutanova, "Bert: Pre-training of deep bidirectional transformers for language understanding," arXiv preprint arXiv:1810.048055, 2018.
[48] L. Schovac and K. Grolinger, "Deep learning for load forecasting: Sequence to sequence recurrent neural networks with attention," leee Access, vol. 8, pp. 36411-364 2020.
[49] S. Holt, A. Hüyük, Z. Qian, H. Sun, and M. van der Schaar, "Neural laplace control for continuous-time delayed systems," " in International Conference on Artificial Intelligence and Statistics. PMLR, 2023, pp. 1747-1778.
[50] M. Schmidt and H. Lipson, "Distilling free-form natural laws from experimental data," science, vol. 324, no. 5923, pp. 81-85 85, 2009.
[51] Z. Qian, K. Kacprzyk, and M. der Schaar, "D-CODE: Discovering closed-form ODEs from observed trajectories," " in International Conference on Learning Representations, 2022. [Online]. Available: http:///eeeeeooooeee////onelvete
[52] K. Kacprzyk, Z. Qian, and M. van der Schaar, "D-cipher: discovery of closed-form partial differential equations." Advances in Neural Information Processing Systems, vol. 36, 2024.
[53] K. Kacprzyk, T. Liu, and M. van der Schaar, "Towards transparent time series forecasting," " in The Twelfth International Conference on Learning Representations, 2024. [Online]. Available: http://peeeeeeeeeeee////oddddd/a
[54] S. Cuomo, V. S. Di Cola, F. Giampaolo, G. Rozza, M. Raissi, and Piccialli, "Scientific informed neural networks: Where we are and what's 's next," Journal of Scientific Computing, vol. 92, no. 3, p. 88, 2022.
machine learning through
[55] S. Greydanus, M. Dzamba, and J. Yosinski, "Hamiltonian neural networks," Advances in neural information processing systems, vol. 32, 2019.
[56] M. Cranmer, S. Greydanus, S. Hoyer. P. Battaglia, D. Spergel, and S. Ho, "Lagrangian neural networks," arXiv preprint wxXiv:2003.06630 2020.
[57] Y. Yin, V. Le Guen, J. Dona, E. de Bézenac, Ayed, N. Thome, and Gallinari, "Augmenting physical models with deep networks for complex dynamics forecasting," Journal of Statistical Mechanics:
Theory and Experiment, vol. 2021, no. 12, p. 124012, 2021.
(58) N. Takeishi and A. Kalousis, "Physics-integrated variational autoencoders for robust and interpretable generative modeling," Advances in Neural Information Processing Systems, vol. 34, pp. 14 809-14 821, 2021.
[59] Z. Qian, W. Zame, L. Fleuren, P. Elbers, and M. van der Schaar, "Integrating expert odes into neural odes: pharmacology and disease progression," Advances in Neural Information Processing Systems, vol. 34, pp. 364-11161111 383, 2021.
[60] A. Wehenkel, J. Behrmann, H. Hsu, G. Sapiro, G. Louppe, and J.-H. Jacobsen, "Robust hybrid learning with expert augmentation." Transactions on Machine Learning Research, 2023. [Online]. Available: http:///eerrcoooooc////////o/cct
[61] C. Geng, H. Paganetti, and C. Grassberger, "Prediction of Treatment Response for Com- bined Chemo- and Radiation Therapy for Non-Small Cell Lung Cancer Patients Using a Bio-Mathematical Model," Scientific Reports, vol. 7, no. 1, p. 13542, Oct. 2017.
[62] I. Bica, A. M. Alaa, J. Jordon, and M. van der Schaar, "Estimating counterfactual treatment outcomes over time through adversarially balanced representations," " in Internationa Conference on Learning Representations, 2020.
[63] N. Seedat, F. Imrie, A. Bellot, Z. Qian, and M. van der Schaar, "Continuous-time modeling of counterfactual outcomes using neural controlled differential equations," arXiv preprint 2rXivv2206.08711 2022.
13

---

[64] V. Melnychuk, D. Frauen, and S. Feuerriegel, "Causal transformer for estimating counterfactual outcomes," " in International Conference on Machine Learning. PMLR, 2022, pp. 15 293- 15 129.
[65] C. C. Kerr, R. M. Stuart, D. Mistry, R. G. Abeysuriya, K. Rosenfeld, G. R. Hart, R. C. Núñez, J. A. Cohen, P. Selvaraj, B. Hagedorn et al., "Covasim: an agent-based model of covid-19 dynamics and interventions," PLOS Computational Biology, vol. 17, no. 7, p. e1009149, 2021.
[66] T. Hiltumen, L. Jones, S. Ellner, and N. G. Hairston Jr, "Temporal dynamics of simple community with intraguild predation: an experimental test," " Ecology, vol. 94, no. 4, pp. 773- 2013.
[67] E. P. Odum and G. W. Barrett, "Fundamentals of ecology," The Journal of Wildlife Management, vol. 36, no. 4, p. 1372, 1972.
[68] V. M. M. Alvarez, R. Roşca, and C. G. Falcujescu, "Dynode: Neural ordinary differential equations for dynamics modeling in continuous control," arXiv preprint arXiv:2009 04278, 2020.
[69] D. E. Rumelhart, G. E. Hinton, and R. J. Williams, "Learning representations by back- propagating errors," nature, vol. 323, no. 6088, pp. 533-536, 1986.
[70] Holt, A. Hüyükk and M. van der Schaar, "Active observing in continuous-time control," Advances in Neural Information Processing Systems, vol. 36, 2024.
[71] S. Holt, Z. Qian, and M. van der Schaar, "Deep generative symbolic regression." in The Eleventh International Conference on Learning Representations, 2023. [Online]. Available:
[72] K. Kacprzyk, S. Holt, J. Berrevoets, Z. Qian, and M. van der Schaar, "ODE discovery for longitudinal heterogeneous treatment effects inference," in The Twelfth International Conference on Learning Representations, 2024. [Online]. Available: http://pppeeeiiinneeeeoouullllpp
[73] T. Liu, N. Astorga, N. Seedat, and M. van der Schaar, "Large language models to enhance bayesian optimization," in The Twelfth International Conference on Learning Representations, 2024. [Online]. Available: https://openreview.nr http://pnneeiieeeee///oummmeeeee
[74] C. Yang, X. Wang, Y. Lu, H. Liu, V. D. Zhou, and X. Chen, "Large language models as optimizers," in The Twelfth International Conference on Learning Representations, 2024. [Online]. Available:
[75] M. Chen, J. Tworek, H. Jun, Q. Yuan, H. P. D. O. Pinto, J. Kaplan, H. Edwards, Y. Burda, N. Joseph, G. Brockman el al., "Evaluating large language models trained on code," arXiv preprint arXiv: 107.:2007.03374, 2021.
[76] S. Holt, M. R. Luyten, and M. der Schaar, "L2MAC: Large language model automatic computer for extensive code generation," " in The Twelfth International Conference on Learning Representations, 2024. [Online]. Available: http:///ppmnewiww.eee/ootntteiee
[77] S. Holt, Z. Qian, T. Liu, J. Weatherall, and M. van der Schaar, "Data-driven discovery of dynamical systems in pharmacology using large language models," in The Thirty-eighth Annual Conference on Neural Information Processing Systems, 2024.
[78] N. Astorga, T. Liu, N. Scedat, and M. van der Schaar, "Partially observable cost-aware active- learning with large language models," in The Thirty-Eighth Annual Conference on Neural Information Processing Systems, 2024.
[79] W. Bonnaffé and T. Coulson, "Fast fitting of neural ordinary differential equations by bayesian neural gradient matching to infer ecological interactions from time-series data," Methods in Ecology and Evolution, vol. 14, no. 6, pp. 1563,, 2023.
[80] S. K. Kumar, "On weight initialization in deep neural networks," " arXiv preprint 17X1v:1004.08863, 2017.
14

---

[81] A. Graves, S. Fernández, and J. Schmidhuber, "Multi MMulti-dimensional recurrent neural networks," International conference on artificial neural networks. Springer, 2007, pp. 549-558.
[82] B. K. Petersen, M. L. Larma, T. N. Mundhenk, C. Santiago, S. K. Kim, and J. T. Kim, "Deep symbolic regression: Recovering mathematical expressions from data via risk-seeking policy gradients," "" in International Conference on Learning Representations, 2020.
[83] S. Hsiang, D. Allen, S. Annan-Phan, K. Bell, I. Bolliger, T. Chong, H. Druckenmiller, L. Huang, A. Hultgren, E. Krasovich et al., "The effect of large-scale anti-contagion policies the covid-19 pandemic." Nature, vol. 584, no. 7820, pp. 262-267, 2020.
[84] O. N. Bjørnstad, K. Shea, M. Krzywinski, and N. Altman, "The seirs model for infectious disease dynamics." Nature methods, vol. 17, no. 6, pp. 202555920000.
[85] Brauer, C. Castillo-Chavez, and Castillo-Chavez, Mathematical models in population biology and epidemiology. Springer, 2012, vol. 2.
15

---

Appendix
Table of Contents
A HDTwinGen Overview
17
B Extended Related Work
17
Benchmark Dataset Enviromment Details
Cancer PKPD
C.2 COVID-19
19
C.3 Plankton Microcosm
19
C.4 Hare-Lynx
20
D Benchmark Method Implementation Details
20
E
HDTwinGen Implementation Details
21
E.1
HDTwin pscudocode
21
E.2
Training HDTwins
21
E.3
HDTwinGen Prompt Templates
22
E.4
HDTwinGen System Requirements Prompts
23
F
Model Optimization Losses
27
Evaluation Metrics
27
H Additional Results
28
H.1
Out-of-distribatioon Experiment and Setup
28
H.2
COVID-19 Unobserved Intervention Experiment and Setup
28
H.3
HDTwinGen top-1 decrcases over time
28
H.4 HDTwinGen Evolution
29
HDTwinGen Ablation Memory
32
H.6 Evaluating Different LLMs
32
H.7 Prompt Ablations with Varying Amounts of Prior Information
32
H.8
Baselines
33
11.9 Procedurally Generated Synthetic Model Benchmark
34
H.10 Interpretability Scale, Performance of only White-Box Models
35
H.1l HDTwinGen Flexibly Integrates Expert-in-the -oopphh-ooopp Feedback
35
H.12 HDTwinGen Accelerates Model Development and Enhances Performance...
36
Hybrid Model Output Examples
38
HDTwinGen can reason about HDTwins
40
16

---

A HDTwinGen Overview
We provide an illustrative example of HDTwinGen working in practice in Figure 5.
User
Generation
Generation 2
Generation
HDTwin
System Rocai 12
Critic L.M
HDTWin
Critic LLM
HDTWin
Actor LLM
Actor LLM
Eyal
Dataset
(Pouis
Actor L.M
Validation MSE
Validation MSF
Validation MSF-
0.0723
Figure 5: HDTwinGen Illustrative Example in Operation. HDTwinGen can generate and further evolve HDTwins for particular system based on user-given system requirements and dataset D Piess of state-action trajectories. First, the system requirements wwhrch include dataset statisties are incorporated into prompt and fed into the modeling agent that returns the code for the HDTwin. This HDTwin is then trained on the training dataset Dimin and a validation loss is computed with Deal- In subsequent generations, the evaluation agent is given the existing generated top- K HDTwins, their corresponding validation losses, and validation losses per component, and is asked to reflect how to improve the HDTwin. This provides detailed, actionable feedback, leveraged from its inherent understanding, and provides this as detailed verbal feedback as H. whereby this feedback is next used with the modeling agent to generate the next HDTwin [P3]. This process iterates several generation times, and the best- rest-performing HDTwin (w.rtt validation performance) is returned. Overall, this produces HDTwin that fulfills [PI-P3].
B Extended Related Work
Sequence models. ML. approaches frequently address system dynamics as a sequential modeling
problem. Initial models like Hidden Markov Models [34] and Kalman filters [35] made simplifying
Markovian and linearity assumptions, later extended to nonlinear settings [36, 37]. Subsequent models, including recurrent neural networks [38], along with their advanced variants [39, 40, 41], introduced the capability to model longer-term dependencies. More recent advancements include attention mechanisms [42] and transformer models [43]. significantly improving the handling of long-term dependencies in sequence data. Another line of work, Neural Ordinary Differential Equations (NODE) [14, 44, 45], interprets neural network operations as integrations of differential equations to model continuous-time processes. Despite being initially driven by natural language processing applications [47], these methods have found utility in modeling complex systems like weather forecasting [46] and energy systems [48]. Furthermore, sequence models can be used in mosel-based RL [70].
Physics-inspired models. Beyond purely data-centric approaches, recent efforts have focused on integrating physical laws into neural system models. Physics-informed neural networks [15, 54] embed physical laws, often as partial differential equations, directly into the learning process. Other notable methods include Hamiltonian Neural Networks [55] and Lagrangian Neural Networks [56], which respect the structural principles of physical systems. These methods are primarily concerned with modeling physics-related phenomenon and require relatively precise knowledge about the system being modeled (e.g. specific differential equations or energy conservation principles) and specialized mechanisms to incorporate them. Regardless, they have demonstrated that the integration of known principles can significantly improve extrapolation abilities beyond the range of training data. We are similarly inspired to incorporate prior knowledge, In contrast, our work aims to integrate more general or partial knowledge flexibly into a hybrid model using LLMs within evolutionary multi-agent framework, while introducing more generalized mechanisms to incorporate loose-form prior knowledge.
17

---

Discovering closed-form models. Closely aligned with our research are techniques aimed at dis- covering closed-form mathematical expressions from data. Symbolic regression [10, 71, 72] and methods like Eureqa [50], SINDY [11], and D-CODE [51, 52] have showcased their prowess in discovering physical laws from experimental observations. However, these techniques can struggle in higher-dimension settings and rely experts to perform the system decomposition to identify the most relevant variables before feeding this information to the algorithm. Additionally, they also rely on experts to specify the function set and mathematical operations that the algorithm uses to search for symbolic expressions. In contrast, our method autonomously learns both the system decomposition and the functional forms of component dynamics, potentially enhancing scalability and efficiency. Moreover, the incorporation of LLMs facilitates the flexible integration prior knowledge at various stages of the search process [73, 74]. Furthermore, using LLMs to generate code, prior work has shown L.MM multi-agent frameworks' ability to excel at large code-generation tasks [75, 76), which could expect to apply here to scale up the size of the generated models in future works. Such future work could also explore acquiring features as well (77, 78].
C Benchmark Dataset Environment Details
In the following, we outline the six real-world system dynamics datasets, where each dataset either a real-world dataset or has been sampled from an accurate simulator designed by human experts.
C.1
Cancer PKPD
Three of our environments that sample dataset from derived from state-of-the-art biomedical hatmmaaakinniic-haaaaadddammmm
Pharmacokinetic- Pharnacodynamic (PKPD) model of lung cancer tumor growth, used to simulate the combined effects of chemotherapy and radiotherapy in lung cancer [61] has been extensively used by other works [62, 63, 64). Here we use this bio-mathematical lung cancer model to create three variations of lung cancer under the effect of no treatments (Lung Cancer), chemotherapy only (Lung Cancer (with Chemo.)), and chemotherapy combined with radiotherapy (Lung Cancer (with Chemo. & Radio.)); for each model we sample respective dataset. First, let detail the general case of Lung Cancer (with Chemo. & Radio.), which comes from the general model (Cancer PKPD Model), and then detail the variations.
Cancer PKPD Model. This is a state-of- ttate-of-the-at -tee-rrtt biomedical Paarmmcokintti-PPPammooyynmmaaa Pharmacodynamic (PKPD) model of tumor growth, that simulates the combined effects of chemotherapy and radiotherapy
in lung cancer [61] (Equation (2))-this has been extensively used by other works [62, 63, 64). Specifically, this models the volume of the tumor srt)) for days f after the cancer diagnosis- where the outcome is one-dimensional. The model has two binary treatments: (1) radiotherapy uj and (2)
chemotherapy ut.
dx(t)) dt
K
(plog
B,C(t)
- -(a,d(t) + B,d(t)²) (,,d(t) )x(()
xx())
(2)
x(t)
Chemotherapy
Radiotherapy
Tumorgrowth
Where the parameters K, p, Ber Oc,,, 9, for each simulated patient are detailed in [61], which also described in Table 3. Additionally, the chemotherapy drug concentration c(I) follows an exponential
Table 3: Cancer PKPD parameter values.
Model
Variable
Paramcter Parameter Value
Growth parameter
7.00 x 10
Tumor growth
Carrying capacity
K
30
Radiotherapy
Radio cell kill (ax
a
0.0398
Radio cell kill (3)
8.
Set s.L. ax ///3= /3=10
Chemotherapy
Chemo cell kill 3e
0.028
decay relationship with a half-life of one day:
de(t) -0.5c(t)
18

---

where the chemotherapy binary action represents increasing the c(t) concentration by 5. 5.0mg/ of Vinblastine given time t. Whereas the radiotherapy concentration d(t) represents 2.0Gy fractions of radiotherapy given at timestep 1, where Gy is the Gray ionizing radiation dose.
Time-dependent confounding. We introduce time-varying confounding into the data generation process. This is accomplished by characterizing the allocation of chemotherapy and radiotherapy as Bernoulli random variables. The associated probabilities, Pe and pr, are determined by the tumor diameter as follows:
(o(1)= = ( ( (D((1-) (0(0)-ō)) 5.))
Pr(() ( (D(0
5,)
D
where Dmax = 13em represents the largest tumor diameter, 0. = 0r = maxx / and (tt) signifies the mean tumor diameter. The parameters To and 7 manage the extent of time-varying confounding. We use 7/-- 7/- = 2.
Sampling datasets. Using the above Cancer PKPD model, we sample N = 1, 000 trajectories, which equates to N = 10, 000 patients, where we sample their initial tumor volumes from a uniform distribution (0) ~ (4 1149). and use the Cancer PKPD Equation (2) along with the action policy of Equation (4) to forward simulate patient trajectories for 60 days, using a Euler stepwise solver. This forms one dataset sample. We repeat this process with independent random seeds to generate Dusin, Dval. Dest Specifically for each benchmark method run for random seed, we re-sample the datasets. For each variation described above, either include the chemotherapy dosing action, chemotherapy and radiotherapy dosing action or neither. We further outline this dataset's system description and variable descriptions with the following prompt template given in Appendix E.4.
C.2 COVID-19
We use the accurate and complex epidemic agent-based simulator of COVASIM [65] to simulate COVID-19 epidemics. This is an advanced simulator that is capable of simulating non-pharmaceutical interventions (such as lockdowns through social distancing, and school closures) and pharmaceutical interventions (such as vaccinations). As this is an agent-based simulator, each agent is an individual in population, and they can be in one of the following states minimally, of being susceptible to COVID-19, 19, exposed, infectious or recovered (which includes deaths). We use the simulator with the default parameters set by the open source implementation of the simulator 5. COVASIM runs simulation for population of individuals. To ensure an accurate simulation, we simulate 24 countries collecting trajectories for each, wherein each simulation we use population size of 000, 000 individuals, and simulate each individual separately (disabling simulation rescaling) and start with a random number of individuals who are infected with COVID-19, T(O) (/(10, 000, 100, 000), and forward simulate the simulation for 60 days. We repeat this process with independent random seeds to generate D Dirain D Diest- Specifically for each benchmark method run for random seed, we re-sample the datasets. We further outline this dataset's '' system description and variable descriptions with the following prompt template as given in Appendix E.4.
C.3 Plankton Microcosm
This describes an ecological model of microcosm of algae, flagellate, and rotifer populations, thus replicating an experimental three-species mey-predator system [66]. We use the dataset made available by [79]S. The dataset consists of a single trajectory of 102 time steps, and we use train, val, test split of 70%. 15% and 15%, ensuring the splits along the time dimension to maintain the integrity of temporal causality, following their chronological order. We further outline this dataset's system description and variable descriptions with the following prompt template as given in Appendix E.4.
'COVASIM is opensource simulator, from which access it here https: com/ stw.uuteeoooisseeeeeooeeeee.ggg covaiii...
'The Plankton Microcosm and Hare-Lynx datasets are both open source and available from https:// github. http:///lleeoooaaeeeeeeeeeeeeeee
19

---

C.4 Hyre-Lynx
This describes real-world dataset of hare and lynx populations, replicating predator-prey dynamics [67]. We use the dataset made available by [79]. The dataset consists of a single trajectory of 92 time steps, and we use train, val, test split of 70%. 15% and 15%, ensuring the splits are along the time dimension to maintain the integrity of temporal causality, following their chronological order. We further outline this dataset's system description and variable descriptions with the following prompt template as given in Appendix E.4.
D Benchmark Method Implementation Details
To be competitive we compare against popular black-box models, which, when modeling the dynamics of system over time, becomes form of ODE model, that neural ODE [14] with action inputs (DyNODE) [68]. We also compare against transparent dynamical equations derived from equation discovery methods for ODEs such as Sparse Identification of Nonlinear Dynamics (SINDy) [11]. Moreover, we compare against the ablations of our method, of the zero-shot generated HDTwin (ZeroShot) and this model with subsequently optimized parameters (ZeroOptim).
DyNODE
DyNODE is black-box neural setwork-based dynamics model [68], that models the underlying dynamics of a system by incorporating control into the standard neural ordinary differential equation framework [14]. We use DyNODE model with 3-layer Multilayer Perceptron (MLP), with hidden dimension of 128 units, with tanh activation functions, and make it competitive by using Xavier weight initialization [80]. To be competitive
the same objective, optimizer and same hyperparameters for the optimizer that we use in HDTwinGen. That of an Adam optimizer [32], with learning rate of 0.01, with batch size of 1,000 and early stopping with patience of 20, and train it for 2,000 epochs to ensure it converges.
Causal Transformer
Causal Transformer is a state-of- thr--rrrttttttttt transformer model for estimating counterfactual outcomes [64]. Due to the complexity of the Causal Transformer, incorporating three separate transformer networks, each one for processing covariates, past treatments, and past outcomes, respectively- which is unique to estimating counterfactual outcomes in treatment effect settings; we implemented only a single transformer model the past outcomes, which is applicable to our datasets and task domains. Specifically, this consists of standard transformer encoder, where the input dataset is normalized to the training dataset. We encode input observed dimension of the state-action into an embedding vector dimension of size 250 through linear layer, followed by the addition of standard positional encoder [64]; this is then fed into a transformer encoder layer, with a head size of 10, dropout 0.1, and the output of this is then fed into linear layer to reconstruct the next step ahead state, of size of the state dimension. We train this model using the AdamW [32] optimizer with a learning rate of 0.00005 and step learning rate scheduler of step size 1.0 and gamma 0.95; also implement gradient clipping to 0.7, with a batch size 1,000 and early stopping with patience of 20, and train it for 2,000 epochs ensure it converges.
RNN
Recurrent Neural Network [81] is standard baseline that is widely used in autoregressive time series next step ahead prediction. We implement this where the input dataset is normalized the training dataset. It consists of a gated recurrent unit RNN taking the state-action dimension in mapping it to a hidden dimension of size 250, with two layers. The output is then fed linear layer convert the hidden dimension back to the state dimension to predict the next step ahead. To be competitive we use the same objective, optimizer and same hyperparameters for the optimizer that we use in HDTwinGen. That of an Adam optimizer [32], with learning rate of 0.01, with a batch size of 1,000 and early stopping with a patience of 20, and train it for 2,000 epochs to ensure it converges.
SINDY
Sparse Identification of Nonlinear Dynamics (SINDy) [11], is data-driven framework that aims to discover the governing dynamical system equations directly from time-series data, discovering a white-box closed-fome foomm mathematical model. The algorithm works by iteratively performing sparse
20

---

regression on library of candidate functions to identify the sparsest most accurate representation of the dynamical system.
In our implementation. we use a polynomial library of order two, which is a feature library of f. = (1, :o; 21, F. }. Finite difference approximations are used to compute time derivatives from the imput time-series data, of order one. Here the alpha parameter is kept constant at 0.5 across all experiments, and the sparsity threshold is 0.02 for all experiments, apart from the COVID-19 dataset where it is set to 1 x 10
APHYNITY
APHYNITY [60] is implemented using domain-specifie expert models as defined in Appendix H.8 combined with a 3-layer MIP, with the same sapeee-prraeterss as in [60].
GP
Genetic programming (GP) is implemented using the implementation and hyper-parameters from the baseline in [82].
HDTwinGen
See the section Appendix E for the implementation details. Specifically, ZeroShot and ZeroOptim ablations of our method using the exact same setup, hyperparameters and prompts. Here ZeroShot generates one HDTwin, and does not fit its parameters, thus evaluating the loss of the model output directly from the LLM. Whereas ZeroOptim, repeats ZeroShot with the additional step of optimizing the parameters of the HDTwin that was generated-agall using the same training as detailed in Appendix E.2.
E
HDTwinGen Implementation Details
Our proposed method follows the framework described in Section 4. We present pseudocode in Appendix E. how the code-generated HDTwins are trained in Appendix E.2, prompt templates in Appendix E.3, system requirements prompts in Appendix E.4 for each dataset, and provide examples of training runs in Appendix J. Specifically, we find a top-K, where K = 16 is sufficient. Additionally, use the LLM of 19T4-1106-Previee with temperature of 0.7.
E.1 HDTwin pseudocode
Algorithm
Pseudocode for Hybrid Digital Twin Generator Framework
Input: modeling context gcontert, training dataset Dreann validation dataset P maximum generations G, top K programs to consider, R 2: Output: Best fitting hybrid model fa.w(0)-.
3: P 0, HI +-0 0
4: for g 1 to G do
5: fa.w(o) LLM H, p(u), Sconteatt) (Generate HDTwin from modeling agent)
6: -())" arg (P) [()e,w()) Drsin) (Fit the model)
7: Compute validation loss per component and overall 6, v
8: p(o+1)= +- pps) @ (Ja, 5,v) (Add HDTwin to the set of top-K HDTwins)
9:
LL (R, pls)) !Generate self-reflection how improve the HDTwin.)
10: end for Return: (o) (The best fitting model that scored the lowest validation loss)
E.2 Training HDTwins
Once the modeling agent has generated an HDTwin fa, it generates it as code. Specifically, it outputs code for a PyTorch [31] neural network module, where this code string is executed, and the module is then trained on the training dataset. The agent importantly observes a code skeleton within its system requirements context Sconteat of which examples of such a skeleton are given in Appendix E.4, However we stipulate that the skeleton must be Modull"" be called
21

---

"StateDifferential", and the parameters must be initialized, and it must define forward function for computing the state differential of the state, where the state and action for that system are input variables to the function. The LIM is instructed to not to modify the code skeleton, only complete it, and return it. This makes it straightforward to process this from text, execute the module, and then train the model.
Specifically,, we train the model on the training dataset, using the standard MSE loss Equation (5), optimizing using the Adam optimizer [32]. We use the same optimizer hyperparameters as the black-box neural network method, that of learning rate of 0.01, with batch size of 1,000 and early stopping with a patience of 20, and train it for 2,000 epochs to ensure it converges, to ensure fair comparison.
Once the model is trained, we compute the val MSE and val MSE per component, which corresponds the val loss per state output dimension Equation (6). Notably, when we append the trained HDTwin back into Pg we include a string representation of it, which includes the values of any named parameters that were initialized in the model. We observe that feeding in the previous optimized parameters helps the LLM in subsequent generations to suggest good starting initial values for the named parameters, from which they can be further refined with the optimization step.
E.3 HDTwinGen Prompt Templates
In the following we detail the prompt templates used. We always use the system prompt when interacting with the LLM.
System prompt
Objective: Write onde to create as effective differential eguatios simulator for given task.
Please note that the code should be fally fanctionsl, No placeholdzrs.
You mast act astonomoddl and you will receive se homn inpai any stage. Yos have reture as axtpst the complete code completing this tadk, and carrectly improve the code to create the most accurate and realistie cimulator possible. You always write est the eode contents. always indent cods with tabs. You cuanot viseclize ay grophical output. exist required.
machine. code cue include
malti-layer perceptions
Use fuactioas proxided. calliag fanctions pooside RFC825) compliast deviation.
without
Defined function schema prompt
"same" "ommllete. StareBifferrnttal code""
"dxerrption"" "Writz out the code body for the Sisteeiifiireeiial tarch model. "
parameters":
Type" " "abject",
preperties I
SSatn oode":
"ype"" string
deseription :. 'Code for the "StateDiffereatial" torch model, inclusive of the model definition, you are zasere.
take your bext guets. must noxempty string.".
I dode_deseription"
"ype"" string".
desceiprion" model."
deseription
Indiesting
only
required t: [' ('SnttDifficreaiial.ddddd dceepiidse
Modeling agent first task prompt
You will system desoription to code differential equation simulator for.
System Description:" [aysteu_ ]escrppccciitionn
Modelling goals:'."
The parameters model optimized observed training dataser the given simulator, genoralize dara.
training daraser very camples,
Requirement Specification:
The code gencraved should achieve lowest possible valication loss of Ie-10 lexs The code generated should interpretable. and fit the dataset accurately possible..
Sceleter cnde to fill Ins [odel
lnow....
You are code crolving machise. and you will called [xenerations) times generate code, improve the achieve lowest possible validation loss.
22

---

differential black
fit perseptrozs). you components
define thrse.
prefarable decompore the gytem into differential equatione (compriimcet) possible.
You
fenctiors. for example Under ao circonstonce on you change the skeletos code function defiaitions. only fill in the code.
The input consors are vectors of shope (batch, size).
Use initially while box models first and then switch to bybrid white and black bos model for the residuals, only after so further best googram iteration improvemeat with while bas modeis.
Make sure your code follows the exect code skeletoa specification.
Use PyTorch.
step-ly stppp generations]
give complete working
genorating Iterationn Jeurrent Iteration)
Reflection prompt
You acncrated following complations. validation you improve
their parameters optimized traising duteact. reflect oa exanples by
your box node) with loss then idd black box models white bot models
programe
( thistory best completions
Here are the top code completions for that you have genereted. sorted for the lowest validation loss lest." (completions)
Please reflect on hos you can improve the code to the dataset as accurately as possible and be iatergretable Think step step-hy- step. Provide only actionable fordhart that has dirzrr changes the code set write our the code, only dexeribe tox inpooved. Where applicable dataset accerately as possible, .900, geserations]
velies optiaiieed parameters reases the cas improved fit
for generating are for nest iteration [iteration] of cself.config
...
Modeling agent in subsequent generations
regenerate the cede function, with the aim improve code to achieve lower validation crrar, the feedback applicable.
penerating you
iteration Igeneration. 10) est Iself contlg generations] iterations something take your paess, generate give
gonerating emply ciring
Please only fill the following code skeleton." ppompte, wwelekke, ceceeiteennnddeeeeeeeeee ende(xeid nmme))
You cannot change the code skeleton. inpet variables..
E.4 HDTwinGen System Requirements Prompts
By following our proposed system requirements format, we constructed prompts for each of the datasets that we evaluated against, which are listed in the following.
Lung Cancer (with Chemo. Radio.)
You will system deseciption to code differential equation simulato: for.
Syseem Description:" Prediction of Treattient Responac for Combined Chens and Radiation Therapy for Son-Small Cell Leag Cancer Paticats Using Bic Mulsematical Model
you mod the cists differential of volums, and chemotterapy drng www.dddgggg_cecdeeeeeeeeeeeeeeee concenttnnn with the inpat actions chenetheropy, dosugs, Seleebeecpp_oosssse
Description
the chemetheragy desage: radiotherapy. docage: Bosage the radiatherapy with snits
chemetherapy
chemoieeraps
aglar3
the chemsiherapy sinblactine with nglne)
The thne wniss is In days.
Additionally these variables have the raspes of: samor_velune: 10.01433, 1176,861] * chemethecapy, dres, concccncnnttiio [0. 9.9975)
chemotheragy chenocherry__ooage dosage: 10,
radisihsrap_ 19, 2.0]
training
potiants,
patiest
days.
Modellisa rouls....
The parameters of the model will be optimized observed training dataxet with the given cimulator.
The oberved training dataact has few amples, and the model mad be ablc zeneralize unsess data.
Requitemntt Spscification=
1 The code geaccated should schieve the lowest possible velidation loss of le or less.
a The code gesecated should le intospretoble. and fi: the dacaset is securately as possible.
Skeietes cods to fill in
class StateDifferentaal Modale):
23

---

Sieeiiiaaaiilcaaiiaal. 1200 Fill the codo
def forward self volumz: torch. Tentor, chemotberagy. conepeooccettaiiie lorch. Texser, cheneiherapy terch Tenoor rod otherpy dosage; tosch.Tessor Tuple[torch Tonsor, torch Teaser]: 1000: Fill is the code here reture voluun dt. L_chemoterrpy.drgggo. www.
Usefal Inow....
You are code cvalving macnize, and you will be called times to generate code, , and improve oode to achieve the lnvest
possible
training componenss (multi-layer patceptrons); towseer, you
salver
patamaters you my differeatial
define
preferable fenetiors, for skelzios ande function definitions. only fill the endz.
(comparineets) possible.
You
Under circomstance en change
The imput tensors vectors shape (batch_size),
Use initially white hox models first and then switch to lybrid white and black hot models for the residusls, only after to further
best program tteration improvement with whirr hax models.
Male sure your code follows the exect code skeletos specification.
Use PyTorch.
Thisk step step-by-step.. then give the complete full working code. generating code for iteration 20:
Lung Cancer (with Chemo.)
You will get system deseription to codz differential equation simulator fer.
System Description:
Prediction of Traammnnt Responee for Combinrd Oene and Radiation Therapy for No-Small Pall ling Cancer Patients Using Mathemattaal Model
you most chencetheragy dosage.
state diffecential of tamor volune, chemotlerray dro. coteenttann;:;;;::; with ingat actions
Description of
valume: Volume the with unlre antt chemotherapy. drwparaappyyyd...oooee conp._.ncceniaaie chemotheragy. dosage:
Concentration the chematherapy drug vinhlaatine with mptmss the chemsiherapy vinblactine with mginr))
The time eniss is in days.
Additionally these variables have the rasges of:
* tamor temor_vplene volone; [0,64195031. 1250:00200509 ctemetherapy. dregyyddre.eeeeeeeeeeeeee concentratiun [0. 9.9915]
chemetheragy 19, 5.0]
training
patiants,
patiest
days.
Modellisa aoals"".
The parameters of the model will be optimized observed training dataxet with the cimulator.
The observed training dataact has few amples, and the model mad be ablc zeneralize unser data.
Requitemntt Specification
i The sode geaccated should achieve the lowest possible validation loss of le less.
a The code acoccated should le intospretable. and fi: the dataset as accurctely as possible.
Skeletes cods to fill
Stareniiiierenntill celf)....imit_.. 1000, Fill the heee
def forward selr
volume:
Tensor, chenetteragy drug coscentration: torch, Teaser, chenetheragy terch. Tensor)
>> Torch torch. Temor. torch. Temoe]:
1000 Fill is the code here
reture (d meme, valume dr. d_chematheraay www.ddunnccncccoeeeeeeeeeo coneecoonttaiiie.....
Usefal to know....
You are code evolving machire, and you will be called times to generate eofe, and impoove the code to achieve the lowest
possible velidation
state differential parametess you my bluck box neural netwerk components (nulii layer perceptross). lowever, you have
solver fil the observed training dataset.
You
Under circametance fenctiors. exampee differential
decompose
(comparimceis)
You
skelsto andr fusctior. defieitions fill the cods
vesiors Use initially white box models first then switch lybrid white and black box models for the residuals, ouly after further
The
(hatch diz).
best program iteration improvemeat with white bex modeis.
Male sure your code follows the exuct code skeletod specification.
Use PyTorch.
Think exep-by-tipp and then give the complete full working code. You are generating code for ineration out 20.
Lung Cancer
You will sysee desoription to coda differentist eguatian cimulator for.
System Description: Prediction of Treatuennt Response for Coubined Cheno and Radiation Therapy for Nor-Small Cell Leug Cancer Patieats Using Bio Mobematical Model
you moit state differential of tanorvelomee. There ase trealments applicd..
24

---

of the
The time enitx it
Additionally these variables have the rasges of:
tenor volune: 10.64195031, 489245534221
The training dataset cossists of 1000 patizats, where each patical is obszeved for 60 days.
Modelling geals:".
The parameters model optimized observed training dataset the given simulator. observed training dateset very samples, the model most
generalize unsees data.
Requrreennt Specification: c..
The code generaved thould achieve lowest possible valication lost of le-d less. The code generated thould interpretable, and fit the dacaset accuratzly possible.
Scaletor cnde to fill is
class StateDifferrnial lon. Module): def i,ii,, sself);
supert super (SteteDifferential
init.__
1000 Fill the here
Tople[torch Tensor
Useful
You code evolving machise, and will be called 20 times to generate code, and improve the code to sehieve the lovest
possible velidation loss.
The model defines the cists differential and will he und with an aot salver to fir the obsrrved training dataser.
You can use my parameters you uses and my black hos neural newerrs components (muli-llayer perceptrons]; nowever, you have to
defins tazse.
preferable to decompose the gysten lnto differential equations coomartiment) possible.
You can use my unary fenctiors for example leg. exp power etc.
Under ao circumstinne cma you change skelztos functios defisitions. , only
cods.
inpul censers vectors shape (batch
Use prugram models
initially
thea improvemest
bybed specification.
box models residuals, oaly ofter
Thiak step by step and then give the complete full working code. You senerating code for iteratinn out of 20.
Hyre-Lynx
You will system descristion code differential equation simulator for.
Syden Description:""
"Modeling Di- Tropropic Prey-Predater Dyzamics is thee and Lyn Ecolugical System
Here you most model the state differential of hare popolatios and lyrs popolation: with the additional lapot of time years. This aime simelate the popalation dysamie within cimplitied di-tropic seotogical syeen compridin pay (hans), and predators (lynxes). interactions include direct predution competition resources. mirroriaz sstural prodator prey mechanisus..
harr_popslation
population:
lyet serving proxy lyns popelation
thousands.
The model should eaptirs the dysamice of thess populations, reflacting the di-trophic pay--prdator interactions and predict popslation sizes based or aistorical dato. The date exhibits 10 yea loss Characteristic oscillatioas dus predator
dynenics.
Additionally these variables have the raspes of: hare popssappaaaioo:: 11.30, 152.65]
lyss population: 1.19, 79.35]
time.in time_in__eee years 11845, 1935]
training dataset coasists time steps validotioz and training dataset coasists time steps each.
Modelling goalertis
The paramciers of the model will be optimized observed training dataxet with the given cimulator, The observed training dataact has few samples, and the motel mad be able generalixe to unses data.
Requirement Specifisation:""
a The code geaccated should sckieve the lowest possible validation loss of le less. a The code aeaccated should le intospretable. and fi: the dataset is accurately is possible.
code
super Diffeenifiiantiall ielt__. 1000. Fill
here
def forward self pare_papulation: torzh, Tensee, lyss_popalation: turch, Temor, tine_it_ycarr: torch. Texsor) Tuple[turch. Tensor. torch. Tenser::
T000 Fill is the code here
rature (4 hare. popalatien dr. d_lyms.popelatinn pwwpppppipo.......
Usefal to know....
You are code crolving machise. and you will called timss to penerate code, and imppove code schieve the lovest possible validation
training
25

---

You
components (auli-llaer perceptrozs): you define decompore the gutrm into differential equationc (compriinett)) poscible.
You prefarable
fenctions, for example skelzios ande fuxction defisitions. only fill the cndz.
Under circumetanee can change
The input censors are vectors of shope (batch,size size).
Use initially white box models first and then switch lybrid white and black box models for the residusls, only after further best poogram iteration snprovemeat with while bas models.
Male sure your code follows the exect code skeleton apecification.
Use PyTorch.
Thisk step sep-ty-step. then give the complete full working
genorating for iteration out 20.
Plankton Microcosm
You will get syete deseription to cods differentisl equation simulator for.
System Desoription:.>
"Modeling Artificial Tri-Trpplcc Pray- Predatar Oscitlations Sinplified Beological Systen
you nost
state diffecential of alpae_pepelation figeellate_eppuaatie and rotifer_ iapet actions This simolate popalation dynamics withia siaplified ti trophic scological system comprisisg plyy (algae). predatios compattioon
predators (flageilaizs), lop predaters (zesferss) mirroring
iningsiid
Description of variables:
population: Total count sige, serving the
pray
peteneddate.pppplatio population: Total count of flagellates. acting internediate predators and prey internediate
pop_predotors, popullaion: Total conat of cotifers. repeesentias predotors
The datasel encapeulates daily population counts seross multiple simulated ceoysiems over period of 100 days, allewing for the analysis of temporal escillations and phase lagi between sprcies.
Additionally these variablex have the rages
prey__epelltiio (0 005095 2.469735]
istermediate, peppllliatt.pppolltton: 00548..3. 1.50000] top predatees population: [0.030316, 0.739244)
The trainigg
validutios add traising coasists
Modelling Josls"".
The parameters of the model will be optimized to an observed training dataset with the given simulator,
+ The oberved traising datssel has very few samples, and the medel mest be able to zeneralize to umsess data.
Requitement Specifization="
The code generated should achieve the lowest postible valication loss, of lesé less code generated should interpretable, and the dataset securately is possible.
Skeletos code fill Module). class StateDifferential def -nint._ self):
cupert Statebifierentiaa ielf.....ttt._..) 1900 Fill the code here
def forward salf pray. populatior: tarch Traser, intzrmadiats, inismadiat.p-eplaaaaaa pepulation: tarch, Temor, top, prodstore, prodatorc.posstttiiee terch. Tenorr Topleftord Tomer Tenshr torch. Teasor]: 1000. Fill is the code here
reture 1d_ pteppey_ppppppeeee_ee___e popslaticn. dt, 4. iiiterncedate, population, dr, ttep__predaters popultttssppsslationn pogulation___d
pocsible avalving
you
genarate impoove
model defines the diffecential will with selver fit the observed dataset.
You
black hot
components (molti-layer peroeptrons): you
deflne these.
It peeferable to decompose the sysen into differeatial equatioss (comparimeets) possible.
You
fanctiots, for eumple lex, exp power etc.
Under an circumetance can you change the skelzios ende fuection defisitions only fill in the code.
The inpat censors are vecior of shape autch_sizei;
Use initially white bus models first and thea switch to bybrid white and black has models for the residusis, only after so further best program iteration inprovemest with white box models.
Make care your code follows the exact orde skeleron specification.
Use PyTorch.
step tep-by-stap.. . give complete working
generating
COVID-19
You will gst system description to code differentiel equation simulato: for.
Systen Description:"
Prediction C0VD-15 Epidemio Dynamics you
differentist enseapitble.
infected interventioss applied. Here
popal_ioo
aosmalized the fixed population.
Dexeription of the variables:
uxscepiile:: Ratio of the pupulatios that is succssiible the virus.
exposed: Ratio of the population that is exposed to the vires,
infectious.
Ixfocred: Ratto of nte popelation that is sctively cartying and traneniting thr simes.
rrcovered: Batic of the population that have recovered from the viras. including those wio ar decrased. popal_oooa Total popoiction of the coustry, , e coastont.
The dine uniss is in days.
Additionally these variables have raspes sxscepitbaa: 1]
26

---

infecsed:
recovered 10.
ppallpppplation [10000, 10008]
The training dataset consists of 24 coustcies, where each country is observed for 60 days.
Modellisz rosis:"..
The parameieis of the model will be optimized to an observed trainias dataset with the ziven simulator.
The observed training dataset has very lew sample,, and the model mast be able genoralixe to unsess dara.
Reqitement Specification
genarared
accurately possible.
Skcletor cndz to fill in:
class StateDifferntial (nn. Modulc):
def init (self):
super (ttarnniiraaaaniiii
calf). Intt............................ O
1000 Fill is the code here
def sell sesceptible: torch, Temor, exgosed: torch. Tonser, infected: torch, Teasor, recovered torch Tensor,
catal_ population: terch,
terch, Tenor)) Tople lorck. Temeer jerch,Texser, lorch, Tensor, lerch. Tensor]
100: the here
sesceprible
exposad
dr)
Useful
You code evolving machite.
will called
and
code schieve lowest
validation loss.
The model defises the differential and will be used with ODE solver to fit the observed training dataset.
You
arameters you vas: defles thise-
black box acural netwerk composents (nulii Eaulii-laye perceptrons): Jowerer, have
ix prefzrable to decompose the sysen into cifferential equations (comparimess) if possible..
You can use my unary fanctiors for esample leg, esp power etc.
Under no oircuisttnce on you change the skeletor code fuectior deffeitions only fill the cods.
The input censors are vectors shope (hatch sixe): box models first
initially
then switch kybrid
bos models residuals, only after
beat poogram iteration improvemest
white box modeis.
sare your follows the exoct code skeletoa specification.
PyTarch.
Thisk crep-hy-step. and
complete full working You genorating for ineration 20.
F Model Optimization Losses
We consider the optimization loss of mean squared error dataset D and also consider higher fidelity mean squared error loss per component.
MSE Loss. Specifically, we optimize the following mean squared error objective,
N T
L(0,D) NXT ΣΣ [1sowwo (2 ((.)))))
(5)
where N x T is the total number of state-action pairs in the dataset. For a given model find the parameters 8" that minimize this loss, i.e. 0r arg ming L(0. Drrai))- Here we optimize & by stochastic gradient descent, using the Adam optimizer [32], however, we note other optimization algorithms could also used such black box optimizers.
MSE Loss per component. We seek to collect detailed quantitative statistics on how well the generated trained system model performs. Therefore, we collect the validation loss per component. Here, we use (f) to indicate the predictions for the jth component.
Noi T
Wm (0" Dvat)
(2) (ts) (t). (t))
(6)
N X
and collect these scalar validation losses per component into vector w (w) w2, Wm]. and compute its mean as v - Ej= W," Doat) i.e. the validation loss.
G Evaluation Metries
We employ mean squared error (MSE) to evaluate the benchmark methods on held-out test dataset of state-action trajectories, denoted as Drest using the loss defined Equation (5) and report this TMSE. Each metrie is averaged over ten runs with different random seeds, and we present these
27

---

averages along with their 95% mean confidence intervals'. For each random seed run, we generate new train, validation and test dataset independently, when we have access simulator. Additionally, when sampling dataset from simulator, we collect datasets of the same number of trajectories for the validation and test set as was used to generate the training set. We then train each baseline on the training dataset and use the validation dataset for early stopping when the method supports this. We then evaluate the performance of each baseline on the test dataset. We repeat this process for each random seed run. We perform all experiments and training using a single Intel Core 19- 19-12 CPU 3, 3.20GHz, 64GB RAM with an Nvidia RTX3090 GPU 24GB.
H
Additional Results
H.1
Out-of-distribution Experiment and Setup
To explore the out-of-distribution shifts, we adapt the Lung Cancer (with Chemo. & Radio.) simulator to have a training dataset of tumor volumes within the test set over all trajectories. To do this we sampled a training dataset, with starting state tumor volumes sampled from x ~ U(,, 574) and then collected 1, 000 trajectories for 60 time steps. We then sampled test dataset out of distribution from that seen in training, with starting initial state of x ~ ((04,, 1149). To ensure that the test state-action states are completely out of the training distribution we slowed down the time interval of the simulator to that of per hour, rather than day, i.e, simulating at At = 24 rather than At = day resolution. We then also verified that the range of cancer volumes seen within the test set do not overlap all with those in the training set, especially throughout and at the end of the trajectory.
range that outside of the tumor volume range
H.2 COVID-19 Unobserved Intervention Experiment and Setup
To assess whether HDTwinGGen can understand and modify its HDTwin, we setup an experiment to see if the modeling agent can adapt already optimized and best-discovered HDTwin for an unobserved intervention, that is a change of the true system state function dynamics, that are latent.
We constructed a scenario with our COVID-19 simulator, to have realistic intervention of lockdown policy, which reduces the physical distance of individuals to one another. In COVASIM we implemented this, where this happens at day 19, dramatically altering how the COVASIM simulator behaves after the intervention is applied. This is known to reduce the effective contact rate B parameter in an SEIR model [83], for the intervention this approximately corresponds to the parameter reducing its value by 75% after the lockdown intervention is applied. Importantly we now sampled training dataset, and validation dataset from the simulator for before the intervention was applied, and only sampled a test dataset after the intervention was applied.
We took the best-found trained HDTwin model for COVID-19, itself created a HDTwin model that an SEIR model with black-box residual components, as seen Appendix I. The agent then was specifically instructed to attempt to adapt this HDTwin model to this unobserved intervention, that was purely described in words. The LLM was able to reason about it, and similarly decreased the HDTwin models internal parameter by approximately 70%, leading to un-observed modelling of such an intervention, when rolling out from the current state, from day 19. This shows the utility of such a hybrid DT framework. Interestingly, the competing methods, such the black-box method, DyNODE and SINDY incorrectly continue the expected trajectory as they are unaware that the underlying system has been intervened on, and its behavior is now different from what they have leart to optimize to.
H.3 HDTwinGen top-1 decreases over time
We observe in Figure 6, that averaged over 10 random seed runs, the top-K, specifically the top-K HDTwin found at each iteration step decreases, implying average the HDTTwin improves in the iteration loop in the beginning stages.
We
use the code at httpss//stackovee com/ httmm//usstioss///e /twoo/nnnnnn//////// from- sample-daaa data compute these
compute- a-
28

---

15
Top
Top 3 Mean
Validation MSE
10
0.0
2.5
7.5
10.0
12.5
15.0
17.5
Number Generations
Figure Average evolution over generations HDTwinGen for the Lung Cancer (with Chemo.
& Radio.)
H.4 HDTwinGen Evolution
Here we list out the specific HDTwins that have been generated at generation steps and are labeled in Figure 3.
Iteration 1. Val Loss: 4,414876937862211 Description: White box only.
Val Loxe: 4.41 [Wherr the val loss dimension volume loss 3.77. chemetherapy concenrration val lock: 6.06)
Iseration:
Ade
Import torch
import torch an as ea
foo typing import Tiple
clase SrateDifferantialinn Medule):
inin sxllf);
super Siateeiifeeretial self)...nil
growth
stghs Paraneter
(0.05:) grouth chemolberap?
(0.03))
def forward tomer, volumz: torch Tentor chemotherapy, wrw.errpydduu_coseetttttttt torch. Teaser, cheneiheragy terch. Tensor,
radiotherapy. dmaze: torch. Tessor) Tuple Temor. turch. Traser]:
Tomor volume change rate rumer volume di 2 self alpha nmar volume self. hera chemotherapy drag concestraation timor volume celf grmm radietherapy. docag. ismor_volume
Chenothrryy drug concestration change rite <_chemetherragt bakeeethrrasy_errgg.oo crup_concntttiiios.... where the dreg administered and then decays
chrmothsrapy dosage ceanmbbraap.duur.ccccccti Assuming simple
d.chematherapy drug.
optimized_ ('olgha" 08451161589117544 "beta" 00233177001549777
00-8021853033003
Iteration Val Loss: 4.233582019805908 Description: White box model with logistic growth for tumor volume and decay rate for chemotherapy drug concentration.
the val per dimessioa volor_voluee
chemetherapy concentration
from typing Teple
class Store Medule):
def in (slff);
self)....niit....
Parameters for the lomer grewtl model
velf alpha m Parameter(ture) densor (0. 13) Timer groutt rate
velf bets m. Parameter( wunch.icneerii.
Tumee death rate chemothersgy
self gamns in.
(torch.
Tomar rate rodiotherapy
kapps 1n.
01) Carrying capacity
envireement chenstherapy dres
Diaa
forwardiself. self.
chenotherapy dugpoosseentttiii: loseh, cheneeterray Temor) Turcc Tentor, turch Teasor]:
radimbeerayy
Logistic modal for voleme
Lumer, volume._dii self-sipha siphs tumor_volame
tamor tamor_velsne volsne xelf, kappa) self, beta
shemotheragy drag concintration tomor volume self gamn radlorherapy dosage tomor volume
Chenotherapy drug conseatration change rale with decsy 6 chheettberaggy wiog_coocoorrrtiio___... chemotheiopy, dosage self.delte chemotherapy drag cone_coonctttationn
return fd_iemor_colmm_____ d_chemotherray www.ooooccnttttiiee__.dd
oprimized parameters l'alpha
'ganna"
"bapp"
29

---

Iteration 3. Val Loss: 2.90204095840454) Description: White box model with logistic growth for tu- mor volume, interaction term for combined therapy effect, and bi-exponential decay for chemotherapy drug concentration.
Val. Loxe: 2.9 (Wherr the val. lost dimeacios ix unservvolmm val loss: 5.36, chemotterapy. wtwperaaappyrreeoooeeeee conceeennnniiticc val loxe: 0.448)
lscration:
impor: torch
inport totch an as 6a
fon typing import Tiple
clase StateDifferential Medule):
soper
growth (0..555335723355553999
Pa.Paareetee wtw.:o......000....
self. bete
xelf Pareaette
sheneiterapy Temor deats dur rdditherapy the eaviromeeni Dreay rate of the chemsiherryy dres calf ara na. Parametar sorch, (essor 013) # Interaction ir briwn chenotherayy and radimbberayy
10 022076775764991183
Traeer death
self Lipps Parameizr (urch.
Parameter torch. (esser(0155.735111211666) (1599735117))))855)
Carrying
xelf Jelia m. Parameter torch densor(.. Senior0.5000000 50000059004444))
def ssel. lomer, tomer_ volums: torch. Temor chemstheergp_dugg dwg.goooooenntiie torch. Teaser, dheaeihergyydossgg terch. Tessor, otheotherypy dosage; torch, Tessor) >> Tuple[torch Tensor torch Teaser]:
Legistie growt model for tomer veleme wilh intoraction tern alpha tumor_volame chemotherapy, drss__cmeenttaiiin chenotherapy drag.
4 Luueer_veluss__
temor_telinc self kappa) self. bela
vomor_ self gams radistherapy__ooga lomor_volume self. ridaahhrryy dasage
decay c. cheeeotheerpp drug
ehematherapy chemotherapy, dosage delta chemotherapy drag. =o
(d temor volume chemotheragy coscentratios
optimized p axcameters l'ulgha" 0.849095224664469 "bela". 0.92880805511566655
'amMIA". 7885955356661 "Kappa"
1032 1578916015625.
delta ": 0. 95336500011199922
"ora ' 0. 700477775S7557551
AdV
Iteration 4: Val loss: 1 1.2533262968663354 Description: White and black box model with logistic growth for tumor volume, Michaelis-Menteee saturation for chemotherapy effect, sigmoid function for radiotherapy effect, resistance development, and black box MLP for residuals.
Val. Loss: 1.25 (Where vol loss dineusioa vonor_volune
loss: 2.49. chemetheropy wrrprddrr_oooooooe concentration val. loss. 0.0203)
ade
import torch
import torch in as na
from syping import Tople
class State Difforenntii (oo. Modole):
def init (self):
super SratrsttatDDffeeeen
self)....itit....
Perameters the lomer grewis model
siphs
grouth cheretherapy
teaser 022840005511444805555
Parameter tensor(1032. (338666111111)))
capacity cavironnent
the chenotherapy drag Interaction between cheneiherapy and radiotherapy self theta- an. Parameter torch demso: (0.53 (0.5:) e Michaclis Mentea satucation ceasiont for chemotherapy
self della
Parameter tensor00, 965661999146992655)
self ua. Parameteri torck. (es.er(00.0024 (00r2..8022188157777888)))
self <ho ta. Parameter; torck. (e.sor(0.5)) () Sigusoid steepness for rediotherapy effect
self 2012 sn. Parsmeter( surch. ()))) # Ressetance development raiz
Black bos compatett for capturing residuats
self recidual alg Soquestialt
na Linearts. 10). ReUU)) (10,
falwardlself, lumes_
toggecsseaattaiie: Teaser, cheneiberapy_oossee terch. Teasor, torch Temsr) Topla Tenoor, torch Trasor]:
volmms with Interactioz iern resieiance Logiatic model for volomm resistance i self.
4 tomor volune__ dt = self dlph tulor_volame i ( tamor tumor self, kappa) (self. thenotherayy dreg. concentration tomor volume) / (self. thete clenoiherapz dreg concentration) - self vel..gamm-
radieiberapy_der
isnor_vsluee 1 ( torch. eattbe.at-eell.dd x radioikerapy_ddoas 133) - self.cta
chemotherupy dresspppy________eeerrrrrsrssss concetaainnn radiotherspy_dougg
tomer_volum: resistance
Bieexponential decay model chemotherapy drog concestration.
6 s.chemotheragy.do wwwgggonncentaaiiin di chrmothrrapy dorage self.deli. chemetherapy wrageppyydyargccccceeooe com.ar.goomcttrrriinnn...
componest
self alspdui__lppierr.. stack (funor(tuur_vveeer chamotberapy coniiccca...tttcccciiiiiiiss chraotkecsy doisge). din=l))
Combite white hot model with rexidusic Lumerr, volume...di residaslel:. 0) 4. Eccherotteerap drug_ cossentrattio. dt residuals [t. I]
rature (d temor volume dr. d chematheragy drig concentration de3
optinized p eacometers l'olpha 0.02999330 'beto". "s 0.27041515707966666,
'2amsea" B0850253743766411. 'kappo"
1032 4199011977775 della". 0.0354118626717755 0.44412924098895. "ett" 0:1205425983335577
0.05195871 458226665
thela 59808 1649681641 "zho
30

---

Iteration 5. Val loss: 0.17562067508897551 Description: White and black box model with dynamic carrying capacity and decay rate, enhanced interaction term, and an expanded black box MLP for residuals.
Val Lover 0.176 (Where the val loss per dimension temor_ val loxs 0.335, chemotherapy. conceetrettttinnnn val loss 0.0163)
Iteration::
from typing Teple
class State Diffeiiffrentittoo Module):
def
init (self):
salf) init 0
Parameters for the lomer grewil model
self alpha in, Paranetzr (torch (tt..//00009599929991010100 000099999399144119922
Tanor arewth caie
telf bets xn. Parameteri tunch (e5.cc000021145551770990000
Tumer death rate due themsiheragy
xelf zamu an. Parameize turch (toooo0...855522377773888888
Tumer desth rate due radieherrpy
self base
Basc carrying capacity of cavironmca:
topps
carrying capacity
7777411866671777997
datay
chemotherapy drig
self Porameter (ercch./2sserr000) Modifier
mod Paromeeter
voluse sheneiherry and radiothecapy chemotherapy
(torch.
6880636666811111)
Michoelis -Menee saturation
self rho sa. Parameter/ torck. Iessert0 44415222560999333
Sigmord steepness for radintherapy effect
xelf 2c1s an. Perametert torch, denser (0.10541796635155))
Resistance development rate
Black hos componert for capturing residusle
self saiidasl_ngg = Seqientialt
na. Liscorte. (4, 20).
na ReE3)))
na Linear(30, 20).
na ReLUU)
na. Lincur (20,
forward self Adjastad eapacity based
satpoitrrrayydrrroossee dro.ercc.ssetttaaiitt Tesser chenetherrpy
Tessor) Tuple[torch
Teaser]:
kappa self. basea_bise self aod (chenotheresy_dougg douge)
Adjasted decay rate based on tamor volume
delta self, deits_bast seif, ded tureersolumm
Logistic grout model for tomer volume with Interactios tern and resistance
rexistance + self, zzta
volums
Lomer volums.... di = relf alpha umar. volume ismor volene lagps) - (azlf. hata chematherapy dww.drruucecteeee cencestration theta chemotheripy wttpbeerrppyuuyy_eeeeeeee (onsssnnsttaaiiin) self-guume rediotherapy, Jasspe tenor_yolous
volese) torch eepcheeeppssseee tho trediotherapy desage 111) self ets chemotthrrpyy cog__coccattattoon radiethcrapy_
model chemstherapy chennttherapy drin
di chrmorhraayy dorage delia tarch marmmmeehhrrayy..o drag concestration, 1.5)
Black box residuel comgoneat
residuals self, residud (eldedllllpptt (tlr(tcc (tack volune. chemotherapy concestiatioon cheaotleruy__ooue radiehhraay dosage). dinal))
Combise white bos model with residuals romer valume... di residasls [ M]
6 .cchemotherag..cog erg._concccntraiinn
di rexiduals (s, I]
reture temor volume__ dt. d.chemotterrayy drug. sogeeeteentiti.o....
optimized_ parametees ('ulpaa 0.2282065226500011, "beta 0.3422962945755. 0.0899328996977651
1031. 203155544444445555 "kappa_ndd" 2.1315995536702.
"deita_base""
0.02670024327
10001002222456-5210555 0.15670098854663)
"eta" 0.00155923882200555.
'cheta" 8.63972267159579
"tho": 535722554550055555
Adv
Iteration 12 Val loss: 0.0722877347745555 Description: White and black box model with logistic growth for tumor volume, tichaelis-Menten saturation for chemotherapy effect, sigmoid function for radiotherapy effect, resistance development, and simplified black box MLP for residuals.
Val
dissension tamor_volene loss: 0.133, chemotherapy wawplayyddddeeeeetttteteeeee concectraation
0.012)
4de
clase State Differettrreeiiialacoa (nc. Medule):
def init (self):
soper (ttatetiifffeeeennn
salf). Inlt 0
Parameters for the lumer growth model
alpha = in Paranetzr(torss (torch ttecettrttccc..eeooo000000 (0.so:(0.03))
V Tamer growth salz
toreeee.ie....e.eeeeeee..eee self bets z.. Parsucterd toreh, (0aser (0.4 4): Tunn death rate due to chemotherapy
telf gemu mn. Parameier (urch (enxor(0000) 0533 Tumer death rate duc in radiothsrapy
self base Parameter teasort-3 01) Base earrying capacity the rasirunment
Parameter/teres (10300)
kappe
drosy carrying eapecity
treatmeent
tessor 13) Madifiar dacay chenatherapy dagg
(0.00)
chemotherapy radiotberapy (10.0/(1000)) Micheelis Mentea saturation
self
chenotherapy
self sa. Parameter jorzh (ensor(0. Sigmosd
radinhhrrapy effect self zeis Parametert turch. teasor (0.15)) Rexistance developmest rate
Black hos componeri for captoring residusts
velf recidual mlg = an. Sequentialt
6a. Lineer (4. 16).
6a LeokkRelUUoo.. 1).
6a Lineert16, 41.
1a LeabyReLUU (0.
1a. Lincur
31

---

forwardiself self,
lesch.Taseer cheaeiberrpy terch.
Tessor) Topla[torch Tensor torch Traser]:
Adjacted apacity based treatment
self kappa_hace xelf mod (chenstheray_ddagg dedaiheray..doug))
Adjasted decay rate based on temor volume
delta szlf, self.deitt__bn bast + deliddela__oo tosee_yolmm:
Logistic growt medel for lomer velime with interastios lera and resistauce resiatance sell. zzia tumer
tumer_colmme colome
Liumer_valene_ a self aipha iumeryvolmme
lumor kagps) (azif hets chemmtherapy, trrpp_.nncsstratioo chenotherrpy coneennnnttttiinn selr.gome rediotherapy dasape tenor, volune tradiotherayy desage 111)
usner volene) (self.the expc-self, rndietheray dasage
chemotherapy drug concestatiion
liucar deeay uodel chemsiherap
t_chenotberrpp wwtgeccceoeeece concentration di chemothrrayy. durage delia turch gon mmnichemthhrrapy.o www. 1.5)
Black bos residasl componest
residaals self residnal nlp (trctt stack ((tumor voleno. chemotherapy drag concestration, chrmotherapy dosage. radietheray desge).. dima)))
Coubite white bos model with residuals
4 Lomor volune di residaels I., M
t_chemoiteraapy wawreetttrryy_cg__ooeeee cog.oooccc........ttriiii....... di residuals [s,
chemotheragy drug
oprimized_
1029. 802490244775 "Koppaod" 2.4545591 286942662,
0.1289799063555005
089927700772244 13.69580192265915.
"rhe 0.2461563389301. 0.1372759033U599077
H.S HDTwinGen Ablation No Memory
Ablation Study, We also ablate HDTwin by removing its memory, only keeping the last hybrid model it generated. We observe decreased performance as shown in Table 4.
Table HDTwinGen Ablation
Lung Cancer
(with Chemo. Radio.)
Method
TMSE 4
HDTwinGen
0.0889:1000453 17.6::2215
HDDTwinGee-nn-mmmoon
H.6 Evaluating Different LLMs
We performed a complete re-run of our main experiments under the same settings, now using different LLM within our HDTwinGen framework, GPT-3.5. These new results are included in Table 5. HDTwinGen can operate fully with a less capable LLM model. However, the generated and discovered models' performance correlates to the underlying L.LM model's performance as expected.
Table 5: Ablation of using different LLMs. Test MSE (TASE) averaged over ten random seeds. HDTwinGen is capable of using other LLM models, however, the best performance results provided with better setter-perfommig LLMs (e.g. GPT-4). The results are presented with + indicating 95% confidence intervals.
Lung
Ling Cherco.) Ling Cancer
Lysx
Macrocom Ture-
COVID-
Medued
(GPPT-4
23r-66
1.22
We also explored the effect of changing the LLMs underlying temperature hyperparameter, here using GPT-4, varying the temperature from 0.7 2.0, where we used 0.7 throughout all our experiments (Appendix E). As outlined in Table 6, we observe that HDTwinGen is still able to operate with different underlying temperature of the LLM.
H.7 Prompt Ablations with Varying Amounts Prior Information
We conducted complete re-run of our main experiments, ablating the prompt, which provides the prior information the form of a textual prior to HDTwinnGen. We provide the ablation results
in Table 7. Specifically the prompt outlined in Appendix E is structured to include separately
32

---

Table 6: Ablation of varying the LLMss temperature. HDTwinGen still operate with different
LLM temperatures. Reporting the test prediction MSE (TMSE) averaged over ten random seeds.
Lung Cancer
Lung Canger (with Chemo.) Lung Cancer (with Chemo, & Radio.)
Hare Lyns
Plankum Sierncosm
COVID 19
Mathod
Yarsg
Tarsz
Taraz
Tarsz 258+:22 29113003 2404-53
Tarnz 4.6e-05-0000000
Yarse 0.0461-..57
HD (Temp 0.0)
8.84:17.8
0.353+-09966 0.085910.0455
0444+11 12
HLTTwinGen (Temp 0.7)
4.4118077 105+2222
0.131 8:43==264
2.1e-06.1.2e 06
1.2212.28
HDTwinGion (Tomp
0,548-600.00
2.32-06433266-000
0.447-3.05
{system description), {skeleton code), and {useful to know) components. Here, prior information is conveyed through the [system description); which describes the system to be modeled, its features, and their ranges; minor system prior information also conveyed through the {skeleton code), as this includes task-specific feature names input variables. The other components of the prompt (e.g. {useful to know)), do not include any task-specific information, and are there to provide general instructions to make the framework work, such as generate pytorch model as code in the response.
We ablate these components of the prompt, first by removing the task-specific prior {system description) labeled HDTwinGen (Partial Context); second, by removing all task- specific priors removing both {system description) and {skeleton code} (where we change the feature names to meaningless names such as x1, x2, etc.) labeled HDTwinGen (No Context); third, by removing only the {useful to know} information that helps the framework, such as instructions to decompose the system, and or combine white-box models with black-box models for the white-box model residuals, labeled HDTwinGen (No Instructions). We observe in the tabulated results (Table 7) that HDTwinGGen can still operate without any task-specific prior information, however having textual prior aids in generating better-performing models, and partially removing HDTWinGGen operation instructions, makes it generate slightly less good models.
Table 7: Prompt Ablations with Varying Amounts of Prior Information. Test MSE TAS averaged over ten random seeds. HDTTwinGen can still operate without any task-specific prior information, however having textual prior aids in generating better-performing models. The results are averaged over ten random seeds, with + indicating 95% confidence intervals.
(WiB Clene.) Lung Cance: (9ith Chans
Lyex
Planin Microcoom
Ture
D550EET33
27715E3
HDTwisGan
30.34:50.3 2,31
2.571-.44
1.52&1.98
2971:519
691e-0644.11000
5.1es100: 1.5ee11
HDltaatinn (No Instrctions)
0212600067
315163
00016000007
17,36472
42113077
0,13IT5,19K
201131C3
X310-60 23:-8
1.21222
H.8
Domain-Specinf Baselines
We performed a complete re-run of our main datasets using domain-specific white-box baselines, as determined by a human expert, as shown in Table Specifically, for COVID-19 modeling, we fit SEIR model [84], Lotka-Volterra model for predator-pres -pray population dynamics (Hare-Lynx & Plankton Microcosm datasets) [85], and a logistic tumor growth model with chemo. radio. effects modeling. HDTwinGen still models the system most accurately, achieving the lowest test prediction MSE on the held-out test dataset of state-action trajectories.
SEIR Model for COVID-19 Modeling. The SEIR model is a compartmental model used in epidemiology to simulate how a disease spreads through population. divides the population into four compartments: susceptible (S), exposed (E), infectious (I), and recovered (R). The transitions between these compartments are governed by the following differential equations:
dS ST -3"
dt
N
dE SI
-EE, dl dl dt dR dt =71,
B N
o F. TI
33

---

Table 8: Table 5. (Test MSE TMSE averaged over random seeds)
Malled Thxees Spare Santean Deorrpraa
Lane OVE Tagal EEOIITTA
Lasor Ial ladie
Plassin Micnema
LOWISS 2018
HiTviatic
44019J0
291-363
132-139
where N is the total population (assumed constant), /3 is the transmission rate, is the rate which exposed individuals become infectious, and is the recovery rate. These parameters are crucial for capturing the dynamics of the disease spread and estimated from data literature.
Lotka-Volterra Model for Predator-Prey Dynamics. The Lotka- Volteraa model describes the dynamics of biological systems in which two species interact, one predator and the other prey. The model is represented by a set of first-order, non- non-linear, differential equations:
dx az -Bxy. dt dy
dxy - dt
where æ and y represent the prey and predator populations, respectively. The parameters a, 8, 7, and ó denote the prey reproduction rate, the predation rate upon the prey, the predator mortality rate, and the rate at which predators increase by consuming prey, respectively.
Lotka-Volte Volterra Triple Species Model. The extended Lotka- Volterra model incorporating a third species involves additional interactions that represent various ecological relationships such as competition, predation, or mutualism. For the sake of illustration, let's consider system with two predators and one prey. The model is described by the following set of differential equations:
dx
dt x(a 3y 8z),
dy y( 1 + cr),
dt
dz
x(-u VI),
di
where: æ represents the prey population. and = represent the two predator populations. o is the natural growth rate of the prey in the absence of predation. 8 and ō the predation rates of the first and second predators on the prey. respectively. 7 and are the natural death rates of the first and second predators, respectively, in the absence of the prey. and are the growth rates of the first and second predators per unit of prey consumed.
Logistic Tumor Growth Model with Treatment Effects. The logistic tumor growth model with chemotherapy and radiotherapy effects incorporates the logistic growth model's capacity to simulate the saturation effect observed in tumor growth, alongside treatment effects. The model can be described as:
dN dt
N K - -
-
r
where N is the tumor cell population, is the intrinsic growth rate of the tumor, and K is the carrying capacity of the environment. C(N) and R(N) represent the effects of chemotherapy and radiotherapy on the tumor cell population, respectively. These treatment functions are often modeled based on dose-response curves and can vary depending on the specific drugs and radiation doses used.
H.9 Procedurally Generated Synthetic Model Benchmark
We performed a complete re-run of our main baselines on a new entirely procedurally generated synthetic model benchmark. Specifically, by procedurally generating synthetic models, this allows us to test how HDTwinGen performs when the LLM has never seen such model. To create diverse synthetic models, we modified the structure of underlying cancer with chemo and radio models to incorporate non-biological random modifications, which include the use of trigonometric operators and division operators. In the following, we provide the exact changes made and the structure the underlying equation.
34

---

Synthetic 1 (inc. sin(w)))) Here the underlying equation is
dr(t)
K ( +, 7isn(wu)) x(1)
ploo
B.C(()
dt
x(1))
Synthetic 2 (inc. -o((t)
dx((t)
(plog
K 3,C(1) (a,,((t) B,,(())) -d((1) x(t)
dl
xx(1)
K Synthetic 3 (inc. log( x)(()))
((1)+N(()))
dx(())
K
plog
x(t) + ((6)) (a, x(t)
Synthetic 4 (inc. ccos(
dx(())
K B.C(1) - (a,,(1)) (, + cos( x(1)
dt p plog
axt))
Synthetic 5 (inc. (C(t)d(t))
dr(t)
K B.C((I - (a, -, x(1)
plog
dl
a(1))
We observe in the tabulated results in Table that HDTwinGen can still generate models that perform well.
Table 9: Procedurally Generated Synthetic Model Benchmark. Test MSE TMSE averaged over ten random seeds. Reporting the test prediction MSE (TuSg) of the produced system models on held-out test datasets across all synthetic datasets. HDTwinGen achieves the lowest test prediction error. The results are averaged over ten random seeds, with indicating 95% confidence intervals.
Synthetic (inc. (emn(ot)
Synihetic (inc. -f(() Syntheti Siine logi soo
Synthetic Line. cooe(
Synthetic (ine #C12581)1
Tuse
Tyrep
Typer
Taras
Twce 13.113.65
DSNODD
6591331
17.6115.3
12.313.32
83:.4,25
SIND,
69:1.87
16:2999722
13.700574
68441143
14.2100 598
ZeneShot
1.18ew0ta22255000
6e90G-24.162-40
4.562-433233 84celB
ZenoOetim
56.2-1
14.5+1117
1 174
56.0++ 12
0.07+1
HDTwinffie
54.2-2.55
0.1707:66.113
1.346-0.377
548--1.9
0.0633-0.1444
H.10 Interpretability Scale, Performance of only White-Box Models
To investigate questions, of how well do the white-box models that HDTwinGen generates perform, we explore an ablation HD'TwinGen where we constrain the generated models to be white-box only, i.e., a mathematical equation with no black-box neural network components. We tabulate this in Table 10, and observe that even when HDTwinGen is constrained to only generate white-box models (HDTwinGen (Only White-Boxy Box)) it still performs acceptably, indicating that the white-box generated models are modeling well the underlying system when fitted to the dataset.
Table 10: Interpretability Scale, Performance of only White-Box Models. Reporting the test prediction MSE (TMSE) the produced system models on held-out test datasets across all benchmark datasets. The results are averaged over ten random seeds, with + indicating 95% confidence intervals.
Lune Canaer
Lang Cancer (wis Chemn)
Lang Cancer (with Cheno. Bualio) Harc-Lanx
Plznktn Micnsosan
COVID-1P
Methed
Zurez
Tarse:
Tarer
Tarez: 33722544
Ture
Toez
HDTwatden OOals White Box)
10.42:00
43
2422200
0000
$0221.7
HDDtweniinn (Whin-los Black-hes meidcals)
K.07
H.11 HDTwinGen Flexibly Integrates Expert-in-the-loop Feedback
Experts play an active role in model development in two main ways:
35

---

1. Initial prompt. Experts can describe the system, and specify modeling goals and task- specific requirements through Soontext (Appendix E).
2. Direct model feedback. In each iteration of HDTwinGen, the expert can provide direct feedback through H to guide model improvement.
To further demonstrate the flexibility of expert involvement in the modeling process, we supply two additional experiments:
1. Expert specifies instructions to produce fully white-box model through Soontext which is provided in Appendix H.10.
2. Expert provides specific feedback on model improvement during the development process through H. We performed this experiment by stopping HDTwinGen during its iterations, providing human expert targeted feedback, and then observing if the subsequent model generated was able to incorporate the feedback by making targeted changes to the underlying model. We confirm that this the case and provide insight that the underlying LLM is able to interpret the HDTwin code model and selectively change parts. We provide figure to illustrate this result, as seen in Figure 7.
P. Expert Focdbbcck
create scaui to preaus charge 2n banor welt guilanty sore citnave aquictess decerperate
Butattrreeiat Soidel
nattethencs effects with neree tors
2,5
Lnsert 30
Uffeer coffuace
Far Juae Lnt
WiTeatten Swns
Figure 7: HDTwinGen can flexibly integrate expert-in-the-loop feedback, if is provided.
H.12 HDTwinGen Accelerates Model Development and Enhances Performance
We seek to determine the runtimes of HDTwinGen model development compared to human-clock time from a human experiment (Mechanical Turk) experiment where participants are asked to refine models and how performant these models are compared to Bayesian optimization (BO) for fixed model budget. To assess this we performed the additional experiments of:
1. Human-driven model development: Hiring two experienced software engineers (Mechani- cal Turks) to develop and refine models. providing them with exactly the same prompts that HDTwinGen uses, using a human instead of the LLM in the model development loop.
2. AutoML: Hyperparameter tuning (HPT) using for DyNODE and SINDy using BO (HPT search space detailed in Table 11).
We performed these experiments on the Lung Cancer (with Chemo. & Radio.) dataset, and the results are tabulated in Table 12, Here we present, the time to generate an individual model (which includes generating the model and training time); the test MSE TMSE after one hour of run-time, and test MSE separately for a budget of only 15 model evaluations.
Analyzing the results provides the following insights:
HDTwinGen Makes Model Development Notably Faster Compared to Human-Clock Time: HDTwinGen takes an average of 45.56 minutes to complete an experiment using 15
36

---

Table 11: HPT details. BO HPT search space for Figure 8.
Method
Hyperparameter Seurch Space
Best Params
DyNode learning_ rate
[Te-5, le-1]
0.0123
DyNode weight_deces
[Ic-5, lc-1]
0.00029
DyNode hidden_ dim
[32. 1024]
788
DyNode aativvvttiv
[tanh, silu, ELU]
tanh normal
DyNode model initialization
[xavier, normal]
SINDy polynomial librarl__rbrrrraae library_degrce degree
[1.3]
SINDy polynomial haammiillllbrr_liiie intraar_interaation
True. False]
True
SINDy threshold
[Ic-5, 1c-1]
0.0194
SINDy alpha
[le-5, Ic-1]
0.0015
Table 12: Method performance comparison. Reporting the time to generate model (in minutes) and the test MSE (TMSE) after one hour of run-time and with budget of only 15 model evaluations. Here indicates 73.1 1% of iteration time is consumed by LLM querying.
Tim tn gereeio mosi [Mrrnee
Tor MSS Turae TAlter mir for iiniinn
Toe MSD Furse: (Riige mnlyy mistel mnstttllaatttttr
mo NONOR
E137
1T12
HPO SiNDy
0.16
13.225
13 23:45
Expers (Mecianicd
HDTwisGon
model evaluations/genecaatinns, whereas human experts took 148.1 minutes (2 hours and 28 minutes) generate and iterate 15 models. Therefore, using HDTwinGen is considerably faster than human-clock time, which is also an advantage of HDTwinGen.
2. HDTwinGen Generates Better Performing Models for a Fixed Model Budget: HDTwin- Gen generates better-performing models compared to the human experts and the Bayesian hyperparameter optimization (HPO) of the baselines of DyNode and SINDy, for a fixed budget of 15 model evaluations. We provide figure showing the best-performing model performance against each generation in a Figure 8. This figure shows that across model evaluations, HDTwinGen still generates models that perform well.
NSDT
HPO for DyNode
102
Human Experts
IIPO for SINDy
Test MSE
101
100
10-1
2
6
8
10
12
14
Number Generations
Figure 8: Best-performing model performance against each generation, for setup in Appendix H.12. HDTwinGen Generates Better Performing Models for Fixed Model Budget.
37

---

I Hybrid Model Output Examples
These are the final discovered hybrid models generated from our method HDTwin, for each respective environment.
Cancer (with Chemo & Radio)
class
StuteDifferetiial (nc. Modulc):
def Jiiit __iiii_...sifi
self) sel()._.iiit
4 Parometere for dee tomor growin andel
self.aipn as Paramcieritores wtwccettttte.teeeeee..ee.. (0.0)) Tanse groech rate
self. beta an. Parameiert
torch (emoor().4)
Touor death zate dus chometheropy
Parametert
rediocherapy
Bore (1C30.))
anstronnznt
celf kappa_mut mod Parametzri
torch.1emoe semor(-220)
Modifier for sarryiag capacity based or treatneui self. delta base
Paramoter
torch. (9.1)
1 Bere Seery rute of dh chunwtherapy draz
mane volume self. delia__o_d nod in Parsneizet
torch. torch-temoon temorr 0003)
Moditier dices
chonenherapy radiotherapy
torch (1D.))
Micharlts -Hnnis
connaut for chematherapy
self.rbo an Parameierr
torch (0.3)
Signald steepners for codtatherapy effect
<elf reta nn. (ircch (omor((0115)) (8.15)) Resietaxes developuent rate
8 Block bus composent for coptering residuals
self. residual mlp (
an Lineari4. 15).
Lr.
def forward sclf.
unor volene: torch, Teasos
chemotherppy drag concestration: torch. Tensor,
chemotherapy. :mage:: torch Temse.
adiotherapy dosage: torch Teaser.
Teplectorch Tessor, torch Tonsor
Adjusted carzing capority based meattver
keppa self, basp
chemstterrpy dosage rediotherepy desoge
Adjusted
based
delta solf,
Lagistic growth mndel for volnus onth Interaction reststonce resistance =1 self.
tenor_voltue
d,tumor volume.
self alpha tumor_volume [1 voisme kapps]
(aelf. beta caemotheraay drog_ coneentnnniiiiiinnnn lomor_volume)
(self. thets
calf gama
radiolterap desage
tenor_voluue
(1
(radiotherapy
addolieeapy
remor_velume
resieiance
Non-liaerr decay model for chrmotherapy drag concentration (_clemotherspyy ciie__cncntaaaiion ciemotheraps_jousg delta terch.pee chcmoliersspy coag_cooscanaatiin 1.5
residaal composerr selt. nlpc
thenotheraps: drsp_ccsccettttiin. chemotherapy dorage. radioihcrapy_ dinal.
Condine white bot model with residvels
Lolmor_ di residuals [t. 0]
Lchemetherspyy breemmthrrspyddre__coo drapspyyddss___eeeeeoee www._ccccttaaiinncccc residualsf:.
38

---

(d_lemon_toiencc__. 4_chzmetherappy www.pppp......ccttttiiiooooo
COVID-19
class StateDafferentialin Medule]:
sppett
paramettoee
(95064197731844
xelf. delts Parameier torch (6.0)3)
Death
# Define xexrel xelmerk for coplarisy congles vallerss Adjust ror srebitecture bazd an grevious Ierttiiii self resddal an n.nn. Sequestiall an Linear(4, 18).
nn ReLl(
nn Lineari12, 135.
an R... Linea(12..
(xelf torch Tessor) Tuple[torch.
total susceplible:
Teasor):
SEI
d L_exposed__.tt seif,bele szscegtible infected self self-sigmo exposed d. Liiffccted_.d1 di self sigme expesed (self gamu self deliaa infected
-celf
eacceptible Infacted
d_rroovered sell gums infecied celf. delts + infected
4 Caleulate randaats anaz the neural sxtwork states torch stack (f enstepiible. exposad, infected, recovered), dima)) residials, cilf rosideal nnt states)
rexidrafs the differsetie! coastioss l_snsceptible__ resiguals 0] d_infecied_ rrcovered di residusie(:.
rntura (d_enseepibbee_.. diespocdd__d d_infested__ d_rreovered__..
Plankton Microcosm
(sclf):
StateDiffeseatial (sit
Define
for the improved model with Prer growth
self self-aiphs Paramzier
(...999911199999555)) Fres deeth rate due predation self pame oa Paramztar[tors (trcc, 1onsar 062177706466730 # tntermealnnee predater offielescy
self. beta Parammiierrtttrrh
self delra na Paramztoriton torch
Top predstor offictency
Curzsiey capacities for end popsiation
self.K peeyyy at. Parcatter(lerd
Carroing cupatity for the pig popslation
self, inteemeeditee n. Porameter leasor (1800005338)))))333335355555
Carrjing capueity for the intermediate prodster
popslatiax
K_top as Paramster htt.///oo..ooor0-888557577555550
Curenisx capecity for top predator populatioe
signa confficieats
Compatiiion carfficisu: for guy pridetors Compatttion corfflcires
predatrr
residvels refiaed orchlrecture
residual
Linear 128:.
an Iropourt 25).
an Linear(12s. 125).
an krL.
an. (b.pout
an Linearc135.
(self.
Tople prey_pogalation: Temor, internediate_popllation Teasor, hog_preditors.psolliii tocch Teasor) olpha grey incernedaate population
population self. self. pray
Intermediste popelation di (xelf
pop_aaopia iniermediate internediate cecaneddee__ppiiiiie popppppoottiio pep_preddtoos__oe pepulaiooo self, signe ptw.aaddiiseppppllioiseoeo
pryypoppsaation celf delta
intormediate popslation self. K. jnitermedditt)
aar_prrcaton__pog popalation d: -self gamu unn_redaters popalation xelf, delts + intermteints www.prddioo__poosliiiii
self.cta tog_prsdator_ populatisse ( top_peedioon_ppsttiiio (elp.K_ip))
Calculate rasidaels aray
residuals ellf renidal (ipttcc poppppptttppppiaiiin intermediate, prarmaiaatt.oppaaaiie. popelaticn, top_predater_.ppplaali dimzl))
Correct differentiel casstioss pray pogelation residvalst..
residsuls
rstura population__.l d_intermediate, L_imeerediatr_ppppittinn popelation__.. d_top d_top_predxttree pwe.atorrcppppiaaiii...
Hare-Lynx
class StateDifferntiill Molole).
def init (self);
super StateDiffrremtial self)... lnit
Define parameters for modified Lotha-Valierea equations
self alpha Paramenetr torch (0..0::0..06)) rate
005))
39

---

delia Pa.aamnei tense:
(0.0021)) Predzter reproductive
Addiiioona porameters rarzying 03) Carrinn copacity bars
aud dar-depprdon offects
eelf
Parameter tenxor(2 tarch.pi 101) Frepsecy orcillations (Mi-grr paried)
xelf
Parameler
self phi ta. toree (e.5or(0.9)) (0.93) Plase skift for exvillations
Neural network for cspturing residnals
self mil = Seqieniialt
i. Li-Lineart: (3, 10).
ia RRUU))
m. Linear(10,
sell population:
1yts popelation: deperdest
years: Tessor) Tuple[torch.
Medified Temse]]
popolattioss
hac_popplltion hase_popslation keppo) self.beta 2ure_popslation
lyss "yssppppalation
populatppittion..oo th -seff
lyss_popsiation xelf.delis her_ppopulaiion lyox_popuiation
iine-dependent excillatory comporent
tine offect a tarch six; selt onep time in yesrs + self selt.phi) phi)
Combie whire box model with ararol aervore rondunte
residauls self residudl__ (ippp (tack (arrppabaepeppllaiie
lyas_populatios, time_eflect), aa=ll)
d_brre_pepellato population dt i- (: d] d_yysx popolction di to residwals [:. 1] return (4 wwwrrrpppslaticc.... dt. __Iro__poplaatio____.1
HDTwinGen can reason about HDTwins
Worked log output of HDTwinGen of part of run, running on the Plankton Microcosm dataset. can reason about structures and parameters.
You will
dexcription code differential equation vimulator for.
Sydee peseriprion:"
"Moddling Artificial Tri-Trphi: Predaarda Oxciliations in Simplified Beolagical Syste
Here you most model the state diffecential of olgae_popelation. flagellate_ population. . and popelt_oo with so iapet actions This sins to simulate the popslation dycamics within simplified ti ni-trophic ecelogical system comprising pray (algae). intermediate predators (fiazeilates). sad les predaturs (rotifer)). The interactions include direct predatios competition resozrses. mitzering asteral intregsild predation nechariums.
Description
prey. pepulaepaiaions: peereeddate.pppolation Total top predators, population: Toual
sigse,
prinary prys
meternediate
and prey
retifers. representing
The datsset encapiulates dails population counts multiple simulated ecosystems period of 100 days, allewing for the anolysis of temporal oscillations and phose lags between sprcies.
Additionally these varisbles have the rage of:
* prey popelatios:: (0 09589, 2.469335]
intermediate popolation: [0, 808438, 1 5000000
* top predators population: 618316, 6.739244
The training dataset coasists time sleps validotios and traising dataset coasists time steps each.
Modelling poals:
The paramzters of the model will be optimized obsrrved training dataser the cimulator, The observed training dataset has very for amples, and the model mod be able generalixe dara.
Requitement Specification:""
a The code gencrated should schieve the lowest possible validation loss of le-6 or less.
a The code geaccated should be interpretable, and fe: the davaset as accurately as possible.
Skeletos code to
class StoteDifferentiall Module):
Fill
def forward self poppppopitt torzh, Temer, Intermediats turch, Temor, iup_predatoor_pppsiaae terch. Tentor)
Toplef iorch. Teneur, Tenur,, torch, Temer]] 1000: Fill is the code here
roture (d pray pogstaticn dr. 4 Intermediate population dt. d top gredators pogstation dt)
Useful to kuse....
You are code evalving machize, and you will be called timss to generate eode, and improve code to achieve the lnvext
pussible validation
state differential my
solver components (nolii-layer perceptrons); towever, you
training
parameters you
peeferable decompose system differeatial equatioss (comparincets)
You
fenctiots
Under circumatance en change skelzios ande fusstion definitions. only fill the ends.
The tensors vectors
(batch
Use initially white bos models first then switch bybrid white and black bus models for the residuals, only after further
best program teeration improvement with whizz bax models.
Makn care your ande fnllows the exact code ckeleron specification.
Use PyTorch.
Think step step-by-step. then give the complete full working code. You generating for iteration 20.
maltiprocessing [User]
40

---

You following
bad
optimized traising examples delineated
improve
validation
Here your previons inerations the best
genarsted. Uee it to if you have echausted whire hot model, L.C. wam white bus model repeate with the ane val loss and then only add black bos models to the white hot modele""
Iteration O, Best Val Loss: 09122467733551
Medel description: white bot nodel
Iteration 1. Best Val Loss: D 0.001299029020011 Medel deseription: white bos nodel with aéditional paremeeorr for compottiion rnd carrying capucity. and ingroved differential equations to better capture the to -iripipiic dynamics.
Here are the top cnde completions to far that you have genersted sorted for the lowast validation lao last:""
Val Loss: 0.0123 (Where the vol per dicension prey poppstoptiiion Toss: 0,0316, intermediate poopulation val loss; 2.13e-05, top predators pepulation lass: 0.00505) Iteration:
torch from typing import Teple
class State Sifeeiifforentiall Module):
del inir (salf):
super StateDifferemiill self)....itii....)
Define the parcaceiois for the Locks Volteria equations
self alpha in Paraneterrr (ood(l)) Prey growth cate
xelf bets z. Parameter( tureh rand Prey death rate to predation
xelf zamu an. Parameter: turch raad
Predater goowh rate due to predatios
Parametzr rand(1)
forward
poyy-popiittio:
poppulationn
Temor, top_predaton puppalatiiip
Tophe
torch. Tensor, torch. Teaserl:
Lolka-Voltert equations
pray wrtpiiite_pppiatioo...
alpht prey pepolation -elf. popalation intermediste popslation sutemecdiate_pppuliiii
self.
prey_population aitemediattspogalaiii self. delta wtp.prettoossppppuaii
d_top- prtp-prddatrs, pepulation_.. -salf prina top- oop-predatere_ popelation celf detra Intermadiate, intermaiaar.ppguiiiii
top_przcaters wow.grrcannr.pppuiana
roture (d proyrropppllattee popslaticn dt. 4. Linterneeiateeppppeeo population,.. di, d_top_predaten pogulation_.... d()
optimized_ excameters ('olpaa". 0.109720100552231 'beta". 0.6870369151474, "smmaa" 0.122403776766653. 'dzlt' 6:821619725113351
Val Loxe 00103 (Whern the val lox dimenstin is pray. poputation val. loxc: 0.00776. intermediate, population val leex: 3.31e 65, top wwt_predators_ppopleee population val. loss: 0,00161) lteration:
impor: torch
import torch an 15 6a
from typing Import Tople
class StoteDifferential Module). (self):
supert StateDifferentiall
Define paramsiers initt.....it.. model
grouth Pray deach dee pridation
velf bots Paramater leasor 10 6870095151172))
kelf
Parameter(tore) tencor(0. (22140557807889333
Fredstor
die to predation
xelf delta
Predutor desth Additional parcineters for competition rnd carrying eapicity
Parametar(turch teneor(0. 6221619725227/58))
self sigme in. Parameter(torch (and(l)) Competition onefficient for pray and intermediate predators
self cia ta. Paraneter; torzk (and
Competition esefficient for intermedizte predators and top predators
self K ns. Parameter ceneeierrtoorhhe (orch.rasdd))...
Cartyiig capacily for the piey popalation
def forward self prey prey.populatios: tarch, Teasee, snrmediattr.ggulattin turch Femor, tup patpprstorr.pppslatiiee terch. Tenor)
Torle Tencer, tarch. Tensor, torch Temer]: prey popelation K) bets pray population
Imporvel differeatial equatioss
pray population poputation
alpha- prey pepotation
Liateceediate__po pttpeceeddieeppppaaioo...
self,
poppppoociiiios intermediaie_ ptt_perdaoossppopplle
self pataeediitteppssatiee popslatios self. delta
ptddp_predator_ppp pepolatttiolli....... jop_prreators poopulatoo self els ww..prdaten_pegsistiie
pama top_peedaters. popelation self. delt, intermediate.
roture (4 pray pogstaticn dr. 4 Intermediste population di. d top predutors pogslation de)
optiized parametere = ('alpha" ts -001S0AASS9851100 "bra'" B6188518600112278 'gmmaas 8 0.56995950011190016 "deln"
0.02177676010233, "signe"; 0,44419259913366 'eta". 044212175545454014 'K'. 0,349550331
er
you can improve
possible interpretable. step -tp--h-- the only deseribe optinized parameters lteration
step- Provide actionall feedback, direct the Do
inpooved. Where applicable the dataset accurately possible. This generating
improved
MuisPrecess! 2024 01 29 06.22:0..022.22 maltigrocessing INFO [Assistant] To insove the cods for iteration 2. consider the followins actionshle Iredback:
1. te velnsorporate Reseures Limitttion for Intermediats Predatore mnd Top Predators so: The currest model inclades a earrying expacity tern (K) for the pray population but does not accosst for resource linitation for intermediate predators and top predators. Introdace carrying ospacity paramsters for these popolations to model their grown saturation.
Refine Interaction Termes; The interaction terms carrestly use simple uultiplicative factors beweera pepulctions. Cossider mice complex interactions that ans isclude poalineas effects or functional cesposses, sech as Rolling Type II fanctional eaptare dynames predution.
resposse
Adjust Negative Gove Rateas: The optimized parameter 'alpha" the proy Molooiiaally unrastistic as it would imply the pray populatton decreasos the abuence predators. Bocure that paramziere contrained positive opuimiration..
negative whick
41

---

Corfficienttse:: Whit compatition canfficlent. Resvalaate the compattion
they hiologically meaniagfel and
the coefficients
during aptiiiatiion.
3. esJatrodcce Dessity Degendent Moctliity for Toy Predstors The model cucrently does not inclede tern for deasity deperddea noriality of top predators. Add lern the differeatial equaiinn for top predators that accounts for intraspecitio competition or other deesity -kgendesi mortality faciors.
Pacontier Initielization Based en Optimized Valses Use the optimized parameter valzes hon ieeratios I as initial valees for the paramzterr steration 1. This on provide better starting paint for oprimiratior process.
Delayses:
ecosytenss cn periods.
between predation changes pepotation charges. improve sigaificantly
gestation
adding drby
below tareshold of mayii-laerr perceptran. the the residuals the white model. This hyrrid approsch potentially model complex interactions that esptured waite box differential equations.
Block Box Modelee; If validotion
bos model improvemente, coseider introducing biack
By implemeating ttese changes. the model choald harom nome hiologically rralietie and my fi the dataser nom accuratlly petcntially cchievias lower validation loss MoiaProcess 2024 01 29 00:02:01.062,
mzltiprocesiing INFO [User]
Please ne regenerate the cede fanctioa. with the aim to improve the cede to ochieve a lower validation error. Use the feedback are penerating
where applicable..
for ileratios something. take your best guess. You have to zenerate code, and caazet give aa cupsyy string aaawer.
total iterations. When gensrating code you unsure
always
following
self): StateDifferentiall ielt_
1000 Fill
def forward self sopulaapoi Tesse,Tessse, istermediate.pspplaaiii Loreh.Tenmee oop.peciaton_pooiiiiie::: terch. Temor) - Toplel Tansor, torch. Tencor, torch Temar] T00E Fill is the code here
reture pospyyyopplitie._... 4_imtermediaie, 4.imemedinn,pppuiatee. popuatiion.. dt, www.:prrcttnn_pssttiiinoo
You czenot change the code skeleton, or lnpet variables.
42

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: The abstract and introduction clearly state the main claims, which are backed up by empirical evidence.
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 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 long is clear that these goals are not attained by the paper.
2. Limitations
Question:
Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification: The main limitations of the work are discussed in Section 8, including the assumption of semantic priors and types of system studied.
Guidelines:
The answer NA means that the paper has 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 (c.g.. independence assumptions, noiseless settings, model sell-speciiiitiinn 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 few datasets or with 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, facial recognition algorithm may perform poorly when image resolution low or images are taken in low lighting. Or speech-to-tex, system might not be used reliably to provide closed captions for online lectures because it fails 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 faimesss.
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 't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community, Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory Assumptions and Proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
43

---

Answer:
[NA]
Justification:
This paper does not introduce any theoretical results.
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 cross- referenced.
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 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.
4. Experimental Result Reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main ex- perimental 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: In Section 6 and Appendices C to G, we thoroughly outline the experimental procedure. This, combined with the released code, will enable all results in the paper to be reproduced.
Guidelines:
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 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 novel architecture, describing the architecture fully might suffice, or if the contribution 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 accomplish this, but reproducibility instructions for how to replicate the results, access to a hosted model (e.g.. in the case of large language model), releasing of a model checkpoint, other means that appropriate to the research performed.
also be provided via detailed
While NeurIPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily a new algorithm, the paper should make it clear how to reproduce that algorithm.
(b) the contribution is primarily new model architecture, the paper should describe the architecture clearly and fully.
(c) If the contribution is a new model (e.g., a large language model), then there should either be way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset or instructions for how to construct the dataset).
(d) We recognize that reprodueibility may be tricky in some cases, in which authors are welcome to describe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should be possible for other researchers have some path reproducing or verifying the results.
5. Open access to data and code
44

---

Does the paper provide open access to the data and code, with sufficient instruc- tions faithfully reproduce the main experimental results, as described in supplemental material?
Question:
Answer: [Yes]
Justification: We provide code at https: btctps://githie com/samhol t/ t/HDTwinGen, Addition- ally, the datasets and preprocessing procedure are described in Appendix C.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurlPS code and data submission guidelines (https://nips.co
cassssonoolli imPolicy) for more details.
While encourage the release of code and data, understand that this might not be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurlPS code and data submission guidelines (https:
/ wttpppccpppp....ppppptpppppppppp wwtpccccppcccceececccccccccctttt cenolic/guidess.
ceabbmmmsssooPPe ssiiisolicyyy for more details.
The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the proposed method and baselines. only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended the paper) is recommended, but including URLs to data and code is permitted.
6. Experimental Setting/Details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: Appendices G details all training and test details, including data splits and hyperparameter tuning procedures for all baselines.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented the core the paper to 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, as supplemental material.
7. Experiment Statistical Significance
Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Question:
Answer: [Yes]
Justification: We report results over 10 seeded runs for all methods, reporting mean and 95% confidence interval of results.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by emor bars, confi- dence intervals, or statistical significance tests, least for the experiments that support the main claims of the paper.
45

---

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.
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 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 range (e.g. negative еπτог 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 text.
8. Experiments Compute Resources
For each experiment, does the paper provide sufficient information the com- puter resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?
Question:
Answer: [Yes]
Justification: In Appendix G, we report all computer resources required to reproduce the results.
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU 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 well 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 into the paper).
9. Code Of Ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics
Answer: [Yes]
Justification: The authors have read the Code of Ethies and confirm that the paper conforms to the code.
Guidelines:
The answer NA means that the authors have reviewed the NeurIPS Code of Ethics.
If the authors answer No, they should explain the special circumstances that require deviation from the Code of Ethics.
The authors should make to preserve anonymity (e.g., s if there is special consid- eration due to laws or regulations in their jurisdiction).
10. Broader Impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [Yes]
Justification: Section discusses the societal impacts of the work.
Guidelines:
The answer NA means that there is no societal impact of the work performed.
46

---

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 direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that 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 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 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 attacks, mechanisms for monitoring misuse, mechanisms to monitor how system learns from feedback over time, improving the efficiency and accessibility of ML).
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification:
This paper does not release any pretrained models or collated datasets that might pose risk to misuse.
Guidelines:
The answer NA means that the paper poses no such risks.
Released models that have 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.
recognize that providing effective safeguards challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Question:
Answer: [Yes]
The paper cites original owners of all code (for baselines) and datasets used in the paper.
Justification:
Guidelines:
The answer NA means that the paper does not use existing assets.
The authors should cite the original paper that produced the code package or dataset.
The authors should state which version of the asset is used and, if possible, include URL.
The name of the license (c.g., CC-BY 4.0) should be included for each asset.
47

---

For scraped data from a particular source (e.g., website). the copyright and terms of service of that souree should be provided.
If assets are released, the license, copyright information, and terms use in the package should be provided. For popular datasets, paperswi thcode. com comdaatasets has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided.
If this information is not available online, the authors are encouraged to reach out to the asset's s creators.
13. New Assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [NA]
Justification:
This paper does not release any new assets.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details of the dataset/code/mmodll as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent obtained from people whose asset is used.
At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsourcing and Research with Human Subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given participants and screenshots, if applicable, as well as details about compensation (if any)?
Answer: [NA]
Justification:
This work does not involve crowdsourcing research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Including this information in the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included in the main paper.
According to the NeurlPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional Review Board (IRB) Approvals or Equivalent for Research with Human Subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
Justification:
This paper does not involve crowdsourcing or research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsoureing nor research with human subjects.
48

---

Depending on the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the 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.
49


Reviews:
[{"review_id": "xIsk6hjXYm", "reviewer": "Reviewer_UNnB", "summary": "An evolutionary search framework for building hybrid dynamics models, especially for so-called digital twins, is proposed. Its notable feature is the use of LLMs for the model proposal and model evaluation in the search. The effectiveness of the method is validated with multiple datasets.", "questions": "As noted in the weaknesses section, the lack of baselines using non-LLM search (e.g., a mere evolutionary search) limits the significance of the experimental results. Could you please elaborate on why the authors did not include such baselines? Or, if the current results already imply something in this direction, it would be helpful if the authors could clarify.", "limitations": "Limitations clearly discussed.", "rating": 6, "confidence": 3, "soundness": 2, "presentation": 3, "contribution": 3, "strengths": "- Learning hybrid digital twins is indeed an important research topic.\n\n- The idea of using LLM for the model proposal and evaluation in architecture search is interesting.\n\n- The experimental analysis is done carefully.", "weaknesses": "Only a major concern is that in the experiment, methods with (evolutionary or whatever) search but without LLMs are not used. As the use of LLMs in the search is the most notable feature of the method, the experiments should analyze such an aspect particularly. The current baseline methods are okay but do not fit this purpose because they are not based on any search algorithms.\n\nA minor thing. Although the paper nicely overviews some of the relevant studies, it seems to lack a series of studies on hybrid modeling in a part of the ML community, for example:\n\n- Yuan Yin, Vincent Le Guen, Jérémie Dona, Emmanuel de Bézenac, Ibrahim Ayed, Nicolas Thome, and Patrick Gallinari. Augmenting physical models with deep networks for complex dynamics forecasting. Journal of Statistical Mechanics: Theory and Experiment, 2021(12):124012, 2021.\n- Naoya Takeishi and Alexandros Kalousis. Physics-integrated variational autoencoders for robust and interpretable generative modeling. Advances in Neural Information Processing Systems 34, pp.14809–14821, 2021.\n- Zhaozhi Qian, William. R. Zame, Lucas. M. Fleuren, Paul Elbers, and Mihaela van der Schaar. Integrating expert ODEs into neural ODEs: Pharmacology and disease progression. Advances in Neural Information Processing Systems 34, pp. 11364–11383, 2021.\n- Antoine Wehenkel, Jens Behrmann, Hsiang Hsu, Guillermo Sapiro, Gilles Louppe, and Jörn-Henrik Jacobsen. Robust hybrid learning with expert augmentation. Transactions on Machine Learning Research, 2023.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "*We thank the reviewer for their thorough review. We are pleased that the reviewer finds our approach interesting, addressing an important research topic, with careful experimental analysis.*\n\n---\n## [P1] Incorporating additional baselines\n\nThank you for this suggestion. While we included SINDy for discovering closed-form governing equations, it doesn't represent an evolutionary search method. To address this, we've added genetic programming (GP) for symbolic regression **[R5]** as a baseline to discover symbolic equations on $\\mathcal{D}\\_{train}$, using implementation and tuned hyperparameters from **[R6]**.\n\nWe have included the additional results in global response **[A1]**, observing that our method discovered superior models. This can be attributed to two key factors: \n1. **Hybrid model design:** HDTwinGen simultaneously optimizes both the specification (functional form) and parameters of mechanistic and neural components. This hybrid approach enables more flexible and powerful modeling compared to symbolic discovery methods, which are limited to purely symbolic equations. \n2. **Efficient search:** By leveraging LLMs, HDTwinGen utilizes domain knowledge, contextual understanding, and learning to enhance search efficiency, leading to more accurate solutions.\n\n**Actions taken:** We have incorporated GP-based SR as an additional baseline in our empirical analysis.\n\n---\n## [P2] Hyperparameter optimization search\n\nTo broaden our comparison with non-LLM methods, we conducted an additional experiment using Bayesian hyperparameter tuning (HPT) search for SINDy and DyNode baselines. For a fair comparison, we matched the number of HPT searches to HDTwinGen's evolutionary search steps (15 precisely). Results presented in **[A3]** demonstrate HDTwinGen's superior performance against these HPT-optimized baselines, highlighting the performance benefits of evolving model specification beyond HPT.\n\n---\n\n## [P3] Extending the literature review\n\nWe appreciate the references to related works, which consider hybrid models combining mechanistic equations with neural components. **[R1]** integrates prior ODE/PDE knowledge into a hybrid model, using specialized regularization to penalize the neural component's information content. An alternative approach investigates performing hybridization in the latent space. **[R2, R3]** consider settings where an expert equation is known, but equation variables are unobservable. Correspondingly, they employ two sets of latent variables: one governed by expert equations and another linked to neural components. **[R2]** introduced specialized regularization to reduce the divergence of the overall model from the mechanistic component, and semantically grounded expert latent variables. **[R4]** leverages the assumptions that expert models remain valid OOD to sample augmented training distributions from expert latent variables. \n\nOur work differs significantly in both motivation and methodology:\n* **Motivational difference:** Our method introduces an automated approach to jointly optimize hybrid model specification (evolving both mechanistic component and neural architecture) *and* its parameters. This contrasts with related works, where experts specify the hybrid model design and optimization is limited to only its parameters. The benefits of our automated approach are substantial: it enables model development at scale, significantly reducing the time and cost associated with human-driven development. Our method involves humans only at key moments (e.g. providing initial context or iterative feedback, if desired), minimizing efforts required from experts. Importantly, our approach has the potential to uncover novel and effective hybrid model designs that might elude human designers.\n* **Methodological novelty:** Our approach uniquely integrates LLMs within an evolutionary framework to discover optimal models, guided by human-provided task context and data-driven feedback. By leveraging LLMs' advanced capabilities in symbolic manipulation, contextual understanding, and learning, we enable the exploration of a vast combinatorial space of hybrid models. This exploration was previously infeasible with standard evolutionary techniques, marking a significant advancement in automated model discovery and optimization.\n\n**Actions taken:** In response to your suggestion, we have *(1)* extended discussions of related works **[R1-R4]** in L269-280. *(2)* We introduced APHYNITY as an additional baseline in global response **[A1]**, where we observed HDTwinGen outperforming APHYNITY. *(3)* We additionally demonstrated HDTwinGen's ability to further evolve a human-specified APHYNITY model in **[A2]**, highlighting the flexibility of our framework in accommodating various degrees of expert involvement through automated optimization.\n\n---\n*We hope that most of the reviewer’s concerns have been addressed and, if so, they would consider updating their score. We’d be happy to engage in further discussions.*", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "QziVM0xW2S", "reviewer": "Reviewer_YMit", "summary": "This paper describes a method to improve digital twins (DTs) with an evolutionary and human driven dynamics, aided by an LLM. The process is meant to optimize hybrid models with both human-driven and computer search-based optimization.The empirical study uses six datasets from the medicine domain. Suitable baselines and ablation tests are also shown. The authors show that a hybrid approach could lead to superior performance with respect to traditional approaches.", "questions": "In section 4, the paragraph Evolutionary Optimization overview: is it possible to provide more details in relation to the computational time and various aspects as per my point 2 above (weakness) in relation to the evolutionary search? \n\nIn the limitations below, I have two questions that can be answered.", "limitations": "The limitations section highlights relying on human input, the knowledge of the LLM,  and the limitation to continuous time systems. However, I feel that the following other limitations should be considered: \n- how can we guarantee that the LLMs, regardless of its knowledge, can provide correct information? Is there a way to test or assess the performance of the LLM? \n- The variability and stochasticity of the evolutionary search could imply that the algorithm might not be consistent across multiple runs. Can the authors speculate on how to measure variability and uncertainty introduced by the evolutionary search?  \n\nBoth points and questions above are particular relevant in critical domains and scenarios such as the medical scenarios used for demonstration in the paper. \n\nI think a discussion on the implications and ethical aspects of using LLMs, human domain knowledge and evolutionary search all in combination for a wide range of domain, e.g. medical domain, is necessary.", "rating": 7, "confidence": 3, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "I have identified the following strengths: \n- The paper addresses an important problem with DT: the need for continuous optimization in front of unseen condition to main the fidelity and effectiveness of DT\n- The method suggests that significant autonomy in the optimization could be achieved by means of an evolutionary process, while still retaining the valuable contribution of human domain knowledge for the setting of the initial parameters and context. \n- The approach is general and could be applied to a variety of problem domains.\n- The approach combines LLMs with optimization techniques such as evolutionary search, producing an advanced algorithm with high potential in the successful use and deployment of DT. \n- The paper provide significant details for reproducibility in the Appendix.", "weaknesses": "1. I don't think that having two agents: a modelling and an evaluation agent makes the system a \"multi-agent\" system. I find it confusing. While this is not technically a weakness, I feel that this system should not be considered a multi-agent system. \n2. As the authors also state evolutionary computation could be computationally expensive and hide several levels of complexity, e.g., how thorough should the fitness evaluation be, how many individuals in a population and how many generations, what kind of mutation is applied and what intensity. I feel that all such aspects are not well described in the main paper despite the evolutionary part of the algorithm is a fundamental part of the algorithm. \n3. Following-up from the previous point, more discussion on the computational complexity of the algorithm in the main paper could be beneficial. \n4. The paper could improve the clarity in relation to which domains, or rather, which domain characteristics are more suitable to the proposed algorithm. It is likely that the algorithm might not perform equally well in different domains\n5. Despite a fair review of existing related methods in section 5, the paper does not state too clearly what are the core novelty aspects that are being introduced with respect to the state of the art: I can infer part of them by carefully examining the cited papers, but it would be better if the authors could highlight the main novelty in a more technical way with specific references to the most similar papers from the literature.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "*We thank the reviewer for their constructive feedback. We’re glad the reviewer finds the approach general and of high potential.*\n\n---\n## [P1] Clarifying \"Multi-agent\" terminology\nWe appreciate your comment on our system's classification as \"multi-agent\". Upon reflection, we agree that our framework doesn't fully embody all characteristics typically associated with LLM multi-agent systems. While our approach incorporates two distinct components (modeling and evaluation) with specialized roles, they are applied sequentially, lacking dynamic, real-time inter-agent communication and autonomous behavior adjustment in multi-agent systems **[R8, R9]**. \n\nWe reclassified our approach as a \"multi-step\" agent, which more accurately reflects its nature:\n* A modeling step that generates and optimizes new model specifications.\n* An evaluation step that assesses generated models reflects on requirements, and provides targeted improvement feedback.\n\n**Actions taken:** We have now updated this in the paper. Thank you for helping us refine terminology and more precisely position our work.\n\n---\n## [P2] Analyzing computational complexity\nThank you for raising this question. Allow us to clarify the evolutionary process and analyze the computational complexity.\n\n**Evolutionary process:** \n* Each iteration proposes and optimizes one new model specification, conditioned on the top-K previous models\n* The population size N=g in the $g^{th}$ iteration\n* Only the new model's fitness is evaluated on the validation set\n* The selection step retains the top-K models\n* Process repeats for G iterations (G=20, K=16 in our experiments)\n\n**Computational complexity.** We introduce lower-case notation to indicate constant time parameters:\n* Evolution/model generation (LLM): $\\mathcal{O}(c)$ (LLM inference)\n* Parameter optimization: $\\mathcal{O}(d)$ (training on $\\mathcal{D}_{train}$, in practice, depends on model/dataset complexity)\n* Fitness evaluation: $\\mathcal{O}(e)$ (on $\\mathcal{D}_{val}$)\n* Model evaluation (LLM): $\\mathcal{O}(f)$ (LLM inference for improvement suggestions)\n* Selection: $\\mathcal{O}(NlogN)$ (top-K)\n\nThe complexity per step is $\\mathcal{O}(c+d+e+f+NlogN)$, with total complexity for G generations: $\\mathcal{O}(G*(c+d+e+f+NlogN))$. We note that constant complexities vary with different datasets and models. For practical insight, we report average wall-clock times on the Lung Cancer dataset: $c=20s$, $d=12s$, $e=0.004s$, $f=20s$. \n\n**Actions taken:** Added complexity analysis and wall-clock times in App H.\n\n---\n## [P3] Discussing domains suitable for the algorithm\nWe appreciate this suggestion and clarify the domain characteristics our algorithm is designed for:\n* Continuous-time dynamical systems with a) Markov property and b) continuous time evolution.\n* Scenarios where hybrid models are likely to outperform purely mechanistic or neural models: a) partial knowledge of underlying dynamics exists, with some non-negligible unclear aspects, b) limited or unevenly distributed empirical data across the state space.\n* Complex model design spaces with multiple hypotheses: where automated exploration can help experts scale model development, review and discover a wider range of potential designs, improving cost- and time-efficiency of model development.\n\nOur experiments focused on systems with partially understood mechanisms and limited observational data. We observed that our algorithm discovered performant hybrid models, leading to more accurate model dynamics, OOD generalization, and improved evolvability with changing conditions.\n\n**Actions taken:** Included a version of this discussion in S8 and in App A.\n\n---\n## [P4] Highlighting technical novelty\nThank you for your insightful comment. Our work introduces the first automated algorithm for designing hybrid digital twins, with the following key innovations:\n* **Automated design:** We present an evolutionary algorithm that optimizes both model specification and parameterization based on an initial expert prompt. This contrasts with existing hybrid models where experts specify closed-form equations and neural network designs, with only parameters being optimized.\n* **Hybrid model discovery:** Our approach leverages LLMs' capabilities in symbolic manipulation, contextual understanding, and learning to explore a vast combinatorial space of hybrid model specifications. This exploration was previously infeasible with standard techniques, as existing discovery approaches are limited to closed-form equations or neural architecture searches within predefined spaces.\n\n**Actions taken:** Added this to the RWs section to better position our work and showcase its novelty.\n\n---\n## Limitations\nWe acknowledge the importance of these discussions. **Actions taken:** We have integrated an extended discussion in App A and a summarized version in S8:\n* **Verification of HDTwin models:** First, HDTwins are composed of human-interpretable components, enabling higher degrees of expert verification. Second, evolved models should undergo rigorous functional testing (e.g., held-out datasets, robustness/fairness metrics) pre-deployment. Future works could incorporate hallucination mitigation strategies (e.g., RAG, constitutional AI, multi-model consensus) to further improve reliability.\n* **Variability in evolutionary search:** One strategy lies in quantifying variability through multiple independent runs and confidence intervals. Another lies in controlling variability by providing tighter requirements/constraints in the initial prompt or iterative expert model steering to narrow the hypothesis space.\n* **Ethical implications:** We recognize the potential for bias propagation from black-box LLMs. To address this, we recommend rigorous review verification for fairness, bias, privacy, and other ethical concerns.\n\n---\n*We hope that most of the reviewer’s concerns have been addressed and, if so, they would consider updating their score. We’d be happy to engage in further discussions.*", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "7IVM2FWx90", "reviewer": "Reviewer_8ECc", "summary": "This paper presents a LLM-powered evolutionary multi-agent algorithm to automate the composition of hybrid digital twins of dynamical systems. Experiments were conducted on several datasets including variants of PKPD model synthesized data, simulated covid-19 data, simulated plankton miscrocosm data, and a real-data of hare and lynx populations. Results versus existing neural ODE, SINDy, RNN, and Transformer demonstrated favorable performance especially in terms of generalizability, sample efficiency, and evolvability of the hybrid model.", "questions": "It’d be helpful for the authors to clarify the relation of the presented work with existing hybrid models, and how it can benefit with the optimization of these hybrid models.\n\nIt’d be helpful if the authors could add details about the search space of the hybrid models, both in the general methodological settings and in each experimental dataset. \n\nFor each experimental dataset considered, it’d be helpful for the authors to add details about the design of the hybrid model (what white box components, what black box components, what search space, to generate what output for what task).\n\n\nClarification on the “spatial” dimension of the experimental datasets would be appreciated.", "limitations": "The authors provided adequate discussion about the limitation and future work for the presented work.", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "Hybrid models are gaining increasing importance for combining white- and black-box modeling. LLM-driven systems for automating hybrid model design, composition, and optimization are critical for enabling the practical adoption of these complex models. The proposed work is thus highly novel and of high potential impact. \n\nThe experimental evaluation considered several interesting datasets and targeted on several important desiderate of hybrid models.", "weaknesses": "While focused on hybrid models, the writing of the paper is missing important related works in recent hybrid models, such as physics-integrated VAE and APHYNITY, in both discussion of related works, and experimentation. \n\nIn experiments, the baselines of the works are primarily black-box models (NODE, RNN, Transformers). As the focus of the paper is the automated optimization of hybrid models, it’d be important to 1) include existing hybrid-models as baselines, and 2) demonstrate the benefits of automated design of these hybrid models (vs. the current approach for optimizing these hybrid models where both the white- and black-box components are predefined and only their parameters optimized). \n\nThe writing of the paper overall is high-level lacking necessary details for properly understanding and assessing the paper, potentially due to the complexity of the method. For instance, methodologically, it was not clear what is the search space of the hybrid model, for either the white-box or the black-box components; experimentally, it was not clear what the hybrid model is trained to do, what is the rough design of the hybrid components, and the MSE metric is being evaluated on what. \n\nThe performance improvements reported (e.g., Table 1) overall seems to be associated with very large standard deviation (in most cases are larger than the mean). Its gain over transformer based approaches is thus not clear given this large fluctuation and the increased complexity of the method.\n\nMost of the experimental data seem to be low in dimension (at each time frame). Please clarify.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "*We thank the reviewer for their thoughtful and helpful review. We’re glad that the reviewer finds our approach to be both highly novel and of high potential impact, though we agree that the addition of existing hybrid models improves the paper significantly.*\n\n---\n## [P1] Extending the literature review\nWe appreciate the references to related works, which consider hybrid models combining mechanistic equations with neural components. **[R1]** integrates prior ODE/PDE knowledge into a hybrid model, using specialized regularization to penalize the neural component's information content. **[R2]** investigates hybridization in the latent space, employing semantically grounded expert latent variables and neural latent variables that are regularized to reduce divergence from the expert component.\n\nOur work differs significantly in both motivation and methodology:\n* **Motivational difference:** Our method introduces an automated framework to jointly optimize hybrid model specification (evolving both mechanistic components and neural architecture) *and* its parameters. This contrasts with related works, where experts specify the hybrid model design and optimization is limited to parameters.\n* **Methodological novelty:** Our approach uniquely integrates LLMs within an evolutionary framework for automated design, guided by expert-provided task context and data-driven feedback. This leverages LLMs' capabilities in symbolic manipulation, contextual understanding, and learning, enabling the exploration of a vast combinatorial space of hybrid models previously infeasible with standard techniques.\n\n**Actions taken:** *(1)* Extended discussions of related works **[R1,R2]** in L269-280. *(2)* Introduced APHYNITY, an existing hybrid model, as an additional baseline in the global response **[A1]**, demonstrating HDTwinGen's superior performance.\n\n---\n\n## [P2] Demonstrating the benefits of automated design of existing hybrid models\nWe agree and have added an experiment demonstrating HDTwinGen's ability to further optimize a human-specified APHYNITY model. Specifically, we seed the process with an expert-designed hybrid model (combining logistic tumor growth model **[R7]** with three-layer MLP). The evolution process is visualized in **[A2]** (in the PDF), showing HDTwinGen further evolving the model specification, resulting in improved performance over the initial model specification.\n\nOur automated approach offers substantial benefits: it enables model development at scale, significantly reducing the time and cost compared to human-driven development, involving humans only at key moments (e.g., providing initial context or iterative feedback if desired). Importantly, our approach has the potential to uncover novel and effective hybrid model designs that might elude human designers.\n\n**Actions taken:** *(1)* Added additional results demonstrating HDTwinGen's ability to further evolve an expert-specified APHYNITY model in **[A2]**, *(2)* Included discussion on benefits of automated design in S8.\n\n---\n\n## [P3] Clarifications on the algorithm\nThank you. We highlight that low-level details are already contained within App E,F,G. Here we clarify each aspect:\n\n* **Search spaces:** We do not impose predefined restrictions on mechanistic models (e.g., primitive/terminal sets in symbolic regression). We also do not specify any architecture primitives (e.g., in neural architecture search) for neural components. The search space is implicitly constrained by the symbolic language, as the evolved model must be valid/executable PyTorch python code.\n* **Hybrid model:** Our hybrid model considers an additive composition of mechanistic and neural components of the form $dx(t)/dt = f(x(t))$, and $f = f_{neural} + f_{mechanistic}$. Evolution optimizes the specification of both mechanistic and neural components, and corresponding parameters.\n* **Training/evaluation and metrics:** The hybrid model is trained jointly to minimize the training MSE loss (Eq 5, App F) using the Adam optimizer. Model evaluation is based on the val MSE and val MSE per component (corresponding to val loss per state dimension, Eq 6).\n\n**Actions taken:** We have highlighted App E more prominently in the main paper.\n\n---\n\n## [P4] Understanding the variability of results\n**Variability.** Thank you for your observation regarding performance variability in Table 1. As HDTwinGen optimizes both model specification and parameters, this is equivalent to searching in a much larger hypothesis space compared to baselines (that only search in the parameter space). The variability stemming from model specifications is also evident in the higher std dev of ZeroOptim results (optimized model based on zero-shot LLM generated specifications). Additionally, despite higher variability, HDTwinGen-discovered models exhibit several key strengths: superior OOD performance (Table 2), sample efficiency (Fig 2) and easier evolvability to changing dynamics (Fig 4).\n\n**Actions taken:** Updated manuscript to discuss variability and the trade-offs of automation explicitly in S8.\n\n**Empirical comparisons with comparable budgets:** We have included additional results comparing performance against DyNODE and SINDY given equivalent budgets. This additional budget is used towards HPT, where each baseline underwent 15 iterations of HPT, matching HDTwinGen's 15 evolutionary iterations. Results provided in **[A3]** (in additional PDF) demonstrate HDTwinGen's superior performance, highlighting the performance benefits of evolving model specification beyond HPT.\n\n**Actions taken:** Included empirical comparisons given comparable computational budget in App H.\n\n---\n\n## [P5] Clarifying spatial dimensions\nThe spatial dimensions of considered benchmarks are Cancer PKPD (4), COVID-19 (4), Plankton-Microcosm (3), Hare-Lynx (2).\n\n---\n*We hope that most of the reviewer’s concerns have been addressed and, if so, they would consider updating their score. We’d be happy to engage in further discussions.*", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "C2imw44gkQ", "reviewer": "Reviewer_jD19", "summary": "The paper presents a neurosymbolic approach to model dynamical systems based on the usage of LLMs and gradient-based optimization. The model perform competitively against the reported baselines. \n\nFirstly, the human modeller is required to define the problem, priors and target metrics in text. Then the first LLM generates python models that are optimized with SGD, subsequently a second LLM distill the results and pass it onto the first LLM which generates a new set of candidate models to be trained again. \n\nPerformance is reported on synthetic data from different models. Proposed model outperform the reported baselines.", "questions": "- Line 55-56: since the work is only shown in the additive case, I suggest you mention this case here also as only showing compositional case results a bit confusing.\n\n- Line 137-139: I find it confusing stating that the outer loss measures generalization. As far as I understand, the outer loss is responsible for model specification, whether this values captures generalization will the depend on whether it's evaluated on the validation dataset, not due to the fact that it is the outer loss.\n\n- I'm unsure whether all the baselines are synthetic data or only some of them. Could you please clarify this?\n\n- In terms of parameter size, how big are the found models $f_{\\theta, w(\\theta)}$ and how big are the baselines models? Both for the algebraic and neural part.\n\n- Line 188: doesn't $P^{(g)}$ contain $w(\\theta)$? If not, could you elaborate on why?\n\n- Why are two LLMs needed? Couldn't a single one do both steps? Have you tried this?\n\n- Suggestion: It'd be interesting to see how the algebraic expression evolve over time during training.", "limitations": "- The model adds a substantial amount of complexity with respect to alternative modelling approaches like the reported baseline transformers. It's unclear that the performance improvements makes it up for the added complexity of the method. \n\n- Approach is evaluated against simulated data. If one of the benefit of the approach is that it can handle small datasets, and generalize well, why not using a real dataset rather than synthetic data? This way benchmarks against SoTA on those benchmark would be standardized and therefor more meaningful that self-reported baselines. E.g. reported results for SINDy are using polynomials of order 2, this choice is not justified in the paper. Or, transformers perform similarly to paper's approach (within standard error), and it is not possible to know whether the choice of training hyper-parameters.\n\n- Like other symbolic approaches, the paper's approach still requires a human to provide modelling priors to the system, hence the claim of line 79 \"a new approach to automatically learn digital twins\" is a bit weaken.\n\n- Besides, Table 10, the explainability potential of the model remains under-investigated.", "rating": 8, "confidence": 4, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "- The model allows to integrate domain knowledge by using LLM\n\n- The proposed model outperform reported baselines.\n\n- Paper is very well written.\n\n- Annex incorporate comprehensive information about the experiments.", "weaknesses": "- It is unclear whether the computational budget of the baselines (in terms of $\\mathcal{L}(f_{\\theta, \\omega}(o), D_{\\text{train}})$ evaluations) is similar to that of the proposed method. \n\n\n- No code is provided.\n\n- The work relies on a closed-sourced model (GPT-4) whose performance also depends on the specific checkpoints used.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "*We appreciate the reviewer’s thorough evaluation and positive feedback.*\n\n---\n## [P1] Clarifying computational budget\nThank you for raising this question. To clarify:\n1. All neural baselines and each HDTwinGen model use identical training settings (2000 epochs, 1000 batch size), consuming the same number of function evaluations per model.\n2. HDTwinGen evolves one new model every generation for 20 generations, using a higher total computational budget.\n3. Exact evaluation counts are provided in the **Questions** response below.\n\n**Empirical comparisons with comparable budgets:** We have included additional results comparing performance against DyNODE and SINDY given equivalent budgets. This additional budget is used towards HPT, where each baseline underwent 15 iterations of HPT, matching HDTwinGen's 15 evolutionary iterations. Results provided in **[A3]** demonstrate HDTwinGen's superior performance, highlighting the performance benefits of evolving model specification beyond HPT.\n\n**Conceptual differences:** HDTwinGen automates specification and parameter design of hybrid digital twins, unlike neural baselines that only optimize expert-specified models' parameters. This automation accelerates and scales model development, improving cost- and time- efficiency. The increased computational complexity trades off with higher automation levels and the potential for better-performing models.\n\n**Actions taken:** We now include the above in App H.\n\n---\n\n## [P2] Code release\nUpon acceptance of the paper, we will release the code, accompanied by extensive reproducibility instructions in App E, to ensure reproducibility.\n\n---\n\n## [P3] Understanding reliance on underlying LLM\n\nWe appreciate this comment and agree that the performance of HDTwinGen depends on the capabilities of the underlying model. To understand how the algorithm scales with different models, we reported the performance of HDTwinGen using LLMs of varying capabilities (App H6). We observed that performance correlated with the capabilities of the underlying LLMs, supporting the hypothesis that HDTwinGen's performance will improve as the underlying models advance. \n\n**Actions taken:** We now include this in S8 of the paper.\n\n---\n\n## Questions\n\n* **Composition operator:** We have revised L55-56 to mention that we focus on additive composition.\n* **Bilevel objective:** Your understanding is correct. The upper-level objective concerns model specification (evaluated on the validation set), while the lower-level objective concerns parameterization performance (evaluated on the training set). This has been clarified in L137-139.\n* **Additional details on benchmark:** The three variants of the Lung Cancer dataset are synthetically generated from PKPD models, while the COVID-19 dataset is generated using an agent-based COVASIM simulator. The Plankton Microcosm (three species) and Hare-Lynx (two species) datasets are real-world ecological datasets. The Plankton dataset is measured in laboratory replication experiments, whereas the Hare-Lynx dataset is measured outdoors.\n* **Parameter count:** On Lung Cancer:\n\n|Baseline|Parameter Count|Seconds per Epoch|# Function (Epoch) Evals|\n|---|---|---|---|\n|DyNODE|33,922|0.02|2,000|\n|SINDy|13|0.01|2,000| \n|RNN|569,002|0.41|2,000|\n|Transformer|2,558,348|0.36|2,000|\n|HDTwinGen|245|0.01|40,000|\n\n* **$\\omega(\\theta)$ in $\\mathcal{P}^{(g)}$:** Your understanding is correct, $\\mathcal{P}^{(g)}$ is the set of optimized models and includes the optimized parameters $\\omega(\\theta)^*$. This has been revised in L188-L190.\n* **Clarification on modeling workflow:** In our framework, the modeling process is divided into two subtasks: model generation (the modeling agent) and model evaluation (the evaluation agent). Each subtask has distinct instructions and performs different roles. Although there are two steps, both are implemented using a single LLM with different prompts and memory for each task. To clarify this, we will describe our system as a 'multi-step' agent rather than a multi-agent workflow. The main distinction is that multi-agent frameworks (e.g., AutoGPT, MetaGPT) involve dynamic interactions and task allocation, whereas our workflow is sequential and fixed. Thank you for this suggestion, which has improved the clarity and presentation of our work.\n* **Evolution of mechanistic component:** We direct you to App H4, where we included the models returned by HDTwinGen in each generation and annotated it with some interpretations of the evolution (including mechanistic and neural components).\n\n---\n\n## Limitations\n\n* **Complexity and performance tradeoff:** We addressed this concern in detail in **[P1]**. As a novel automated hybrid design strategy, HDTwinGen's computational complexity is indeed higher than manual model building, similar to other automated approaches (eg AutoML, NAS). However, this complexity is justified by: improved model performance, and enhanced scalability in model development, potentially reducing overall cost and time in modeling. This is further supported by our new results (**[A3]**), highlighting that given comparable budgets, HDTwinGen discovers better-performing models.\n* **Standardized comparison:** Please see our response above clarifying synthetic vs real-world datasets. Regarding your comment about the performance benefits of HDTwinGen vs HPT of neural baselines, we addressed this with additional empirical results with comparable computational budgets, which we discussed in **[P1]**.\n* **Sharpening claims:** We have refined our claim in L79 to \"a new approach to automatically design digital twins given human modeling priors.\"\n* **Explainability:** While we have provided some preliminary interpretations of discovered models in App H4 and App H10, we agree that this is an important future direction. We have highlighted this by revising Section 8.\n\n---\n*We hope the reviewer’s concerns are addressed and they will consider updating their score. We welcome further discussion.*", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:
*We are grateful to the reviewers for their insightful feedback and constructive comments that have improved the paper.*

We are encouraged by the reviewers' recognition of our work's novelty and potential impact. Reviewers highlighted our approach as "highly novel and of high potential impact" (**8ECc**) with an approach that "integrate domain knowledge by using LLM" (**jD19**) and achieves "significant autonomy in the optimization by means of an evolutionary process" (**YMit**). They noted HDTwinGen as an "advanced algorithm with high potential" that is "general and could be applied to a variety of domains" (**YMit**).

We are pleased that the reviewers agreed that our work addresses "an important research topic" (**UNnb**, **YMit**) and "LLM-driven systems for automating hybrid model design, composition, and optimization are critical" (**8ECc**). Regarding our empirical analysis, reviewers commented, "the analysis is done carefully" (**UNnb**), "outperforming reported baselines" (**jD19**), and "targeted several important desiderata of hybrid models" (**8ECc**).

We have also taken the reviewers’ feedback into account and made the following key changes to improve the paper:

* **[A1] Additional baselines:** We have added comparisons against *(1)* APHYNITY and *(2)* Genetic Programming for symbolic regression.
* **[A2] HDTwinGen optimization of an existing model:** We have provided an additional experiment demonstrating HDTwinGen's ability to further optimize human-specified APHYNITY models used to seed the evolutionary optimization. This highlights the flexibility of our method in accommodating various degrees of expert involvement and its ability to automate the optimization of an existing hybrid modeling technique.
* **[A3] Empirical results with comparable computational budget:** We further compared HDTwinGen against *(1)* DYNODE and *(2)* SINDy with comparable computational resources, allowing 15 iterations of hyperparameter optimization for baselines and 15 iterations of evolutionary model optimization for HDTwinGen.
* **[A4] Extended related works:** Following concrete recommendations, we have expanded the related works section to discuss and compare against existing hybrid models, including Physics-VAE, APHYNITY, LHM, and AHMs.

These revisions have been reflected in the updated manuscript, with additional empirical results provided in the attached PDF.

We believe these updates and our individual responses address the reviewers' concerns and strengthen our paper. We remain open to further feedback.

With thanks,

The Authors of #17161

---

### Additional References

**[R1]** Yin, Y., et al. Augmenting physical models with deep networks for complex dynamics forecasting (2021)

**[R2]** Takeishi, N. and Kalousis, A. Physics-integrated variational autoencoders for robust and interpretable generative modeling, (2021)

**[R3]** Qian, Z., et al. Integrating expert ODEs into neural ODEs: pharmacology and disease progression (2021)

**[R4]** Wehenkel, A., et al. Robust hybrid learning with expert augmentation (2022)

**[R5]** Koza, J.R., Genetic programming as a means for programming computers by natural selection (1994)

**[R6]** Petersen, Brenden K., et al. Deep symbolic regression: Recovering mathematical expressions from data via risk-seeking policy gradients (2020)

**[R7]** Vaghi, C., et al. Population modeling of tumor growth curves and the reduced Gompertz model improve prediction of the age of experimental tumors (2020)

**[R8]** Hong, S., et al. MetaGPT: Meta Programming for A Multi-Agent Collaborative Framework (2024)

**[R9]** Significant-Gravitas. (2023). AutoGPT. GitHub repository. https://github.com/Significant-Gravitas/AutoGPT

Meta-review: The paper proposes a neurosymbolic approach to modelling complex dynamical systems, also using domain knowledge and LLMs. 

All reviewers recommended acceptance. Reviewers appreciated the novelty of the approach. The authors added baselines in response to reviewers' concerns.


Therefore, I recommend acceptance.
Decision: Accept (spotlight)