Automatically Learning Hybrid Digital Twins of Dynamical Systems

Authors: samuel  holt
tennison  liu
mihaela  van der schaar
modeling  context
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-15 14:54:26

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.

Main paper: 
 introduction: 
 Digital Twins (DTs) are computational models that accurately simulate the states and temporal dynamics of real-world systems [1, 2] . They are 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 a contagious disease containing various components, including infection rates, recovery rates, population movement, and intervention strategies. DTs 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. to curb disease transmission) [5, 6] .Desiderata. A notable differentiator between DTs and general machine learning (ML) models is the emphasis on generalization. DTs are designed to simulate completely unseen scenarios or interventions at inference time. Therefore, a crucial consideration is ▶ [P1] out-of-distribution generalization: the ability to generalize to state-action distributions beyond those observed during training. This challenge is often compounded by the scarcity of observational data available to accurately 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 epidemiological modeling and treatment planning, where DTs are regularly updated to reflect fundamental changes in disease 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 DTs primarily utilize two approaches: mechanistic models or MLbased neural models. Mechanistic models, denoted as f mech , 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, f neural , 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-action distributions, but struggle in 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 Digtal Twins (HDTwins) combine the strengths of both approaches through compositions of neural and mechanistic components, i.e. f = f mech • f neural . Here, f mech 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, f neural 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 problem-specific models [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, 21] . To address this, we introduce HDTwinGen, a 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, HDTwinGen efficiently evolves better performing hybrid models with informed modifications.Contributions: 1 Conceptually, we present the first work in automated hybrid model design, jointly optimizing model specification and parameterization of hybrid digital twins. 2 Technically, we introduce HDTwinGen, 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 DTs, achieving ▶ better out-of-distribution generalization, ▶ sample-efficient learning, and ▶ increased flexibility for modular evolvability. 
digital twins of dynamical systems: 
 A dynamical system S := (X , U, Φ) is a tuple of its d X -dimensional state space X ⊆ R d X , an (optional) d U -dimensional action space U ⊆ R d U , and a dynamics model Φ. The state at time t ∈ T ⊆ R + is represented as a vector, x(t) ∈ X and similarly the action taken is represented as a vector u(t) ∈ U. The continuous-time dynamics of the system can be described by dx(t)/dt = Φ(x(t), u(t), t), where Φ : X × U × T → X . We optionally consider the existence of some policy π : X → P (U) that acts on the system by mapping a state x(t) to a distribution over actions u(t). Digital Twins. Digital twins (DTs) aim to approximate Φ : X × U × T → X using a computational model f θ,ω(θ) ∈ F learned from data. Here, we use θ ∈ Θ to denote the specification of the model (e.g. linear) and ω(θ) ∈ Ω(θ) to indicate the set of parameters specified by θ. Additionally, F, Θ, and Ω(θ) are the spaces of all possible models, specifications, and parameters, respectively. Next, we outline the key desiderata for a DT:[P1] Generalization to unseen state-action distributions. As DTs are required to simulate varying conditions, they should extrapolate to state-action distributions not observed during training time. Formally, the generalization error E (x(t),u(t),y(t))∼p OOD [L(f θ,ω(θ) (x(t), u(t)), y(t))] should be minimized, where L is some loss function, and p OOD represents the out-of-distribution scenario.[P2] Sample-efficient learning. Given the often limited availability of real-world data, DTs 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. θ and ω(θ) should be easily adjustable to reflect changing system dynamics. For the purpose of model learning, we assume access to an offline dataset containing N ∈ N + trajectories, where the measurements of the systems are made at discrete time points[T ] = [t 1 , t 2 , . . . T ]. This dataset, D = {{(x (n) (t), u (n) (t), y (n) (t)) | t ∈ [T ]}} N n=1, contains state-action trajectories sampled regularly over time, where y (n) (t) = x (n) (t + ∆t) represents the subsequent state. 
hybrid digital twins: 
 HDTwin. A Hybrid Digital Twin is a composition of mechanistic and neural components, represented as f θ,ω(θ) = f mech • f neural [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 to adapt efficiently to remain accurate in changing conditions ([P3] ). In this work, we focus on additive compositions, f θ,ω(θ) = f mech + f neural , as they are more interpretable. Additionally, it enables individual contributions of mechanistic and neural components to be 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 HDTwins [29] .Learning the hybrid model can be decomposed into two steps: (1) model specification, or learning the structure, θ ∈ Θ, of the dynamics function that describes how the system evolves over time; and (2) model parameterization, which estimates the specific values of parameters ω(θ) ∈ Ω(θ) for a given specification θ. For instance, the logistic-growth 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: L inner (θ, ω(θ))Here, the upper-level problem involves finding the optimal specification θ * that minimizes the outer objective L outer , while the lower-level problem involves finding the optimal parameters ω * (θ) for a given specification θ that minimizes the inner objective function L inner . To be more concrete, the M U K I S l 6 U V U K R K X K T X J 6 1 / 8 U N Y J 4 3 + j o t S T A q Y a Z l L D u i l a f c n K w D n H F Q 9 b K 7 q / i x q 6 G f K 6 r 1 + P q 0 Z z g X C P m W m E D P o r 2 h 0 9 b Z N T K L G G 5 l Q C H 9 Y V T q p j F 7 x a J + x z h 5 T m U F H W / i / j m e 3 O p 5 t d P Q 8 Y s 2 0 2 4 s H 8 T L o X Z C s Q Y + s Y z j t 3 r D M 8 K o Q G r k C 5 8 Z J X O K k B o u S K 9 F 0 W O V E C f w a Z m L s o Y Z C u E m 9 H G t D 3 3 g l o 7 m x / m i k S / X v i h o K 5 x Z F 6 j P b I b p N r x X / 5 Y 0 r z D 9 O a q n L C o X m q 4 v y S l E 0 t N 0 R z a Q V H N X C A + B W +EQUATION✓,!(✓)⇤ , (1) , (1) ), . . . 3) The HDTwin is evaluated based on model loss υ and component-wise loss δ. Subsequently, the model pool P (g) is updated with top-K models. 4) The evaluation agent provides targeted feedback for model improvement H (g) by analyzing models in P (g) using performance metrics requirements outlined in S context . This iterative loop repeats for G iterations.(f (K) ✓,!(✓)⇤ , (K) , (K) )}outer objective measures the generalization performance, empirically measured on the validation set L val , while the inner objective measures the fitting error, as evaluated on the training set L train . Combinatorial search space. The space of possible specifications Θ (e.g. different networks, functional forms) is discrete and combinatorially large, while Ω(θ) represents the continuous space of parameters to be optimized. Selecting the optimal θ, ω(θ) 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, our work is the first to address the problem of automatic HDTwin development, where we incorporate LLMs (combined with offline optimization tools) to automatically optimize both the specification and the parameterization of hybrid models. 
hdtwingen: automatic design of hdtwins: 
 Human experts craft models by making strategic design decisions based on their domain knowledge, starting with a 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 domain-consistent knowledge into the model specification, particularly through the symbolic representation f mech ; ▶ efficient search, by making intelligent modifications to the specification to converge more efficiently on the optimal hypothesis; and ▶ contextual understanding, enabling the algorithm to incorporate task-specific context and targeted feedback for model improvement [22, 23, 24] .Overview. We operationalize this insight through HDTwinGen, an evolutionary algorithm that iteratively evolves a population of candidate solutions to automatically search for the best HDTwin. Our approach employs a framework comprising three key elements: (1) human experts provide an initial system description, modeling objectives, and requirements as a structured prompt; (2) a modeling agent proposes new model specifications, optimizes their parameters on a training dataset, and collects validation performance metrics; (3) an evaluation agent assesses the proposed models using both data-driven performance metrics 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.1.The optimization process begins with a human expert providing a structured prompt, referred to as the modeling context S context . 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. 2. The modeling objective specifies quantitative performance requirements via a metric L. 3. 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 mechanistic 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, S context includes a skeleton code to guide the synthesis of executable code in a predetermined format. For illustrative purposes, an example of S context is provided in Appendix E.4.Evolutionary optimization overview. Given S context as input, HDTwinGen performs G iterations of optimization, where G ∈ N + . The population of proposed HDTwins at iteration g is represented as P (g) . Each iteration creates a new candidate model based on previously created models in P (g) 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 g∈[G] |P (g) | = K. Each model in P (g) 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 (G) is selected as the final model.Proposing HDTwins. 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 ▶ P (g-1) : the set of top-K previously generated models; ▶ H (g-1) : the most recent feedback produced by the evaluation agent (where on the initial step, g = 1, both are empty, i.e., H (0) = ∅, P (0) = ∅); and ▶ S context : the modeling context. The modeling agent generates a model specification θ 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 of in context-learning and evolve improved specifications in subsequent generations [22] . Each generated specification emits its corresponding parameters, ω(θ) are fitted to the training set D train . More formally, we represent this generative procedure as f θ,ω(θ) ∼ LLM model (H, P (g) , S context ).Model specification. To generate model specifications, the modeling agent decomposes the system into a set of components, with each component describing the dynamics of a specific state variable. In other words, for a system with d X state variables, there will be d X components. Each component is characterized by its own set of inputs and a unique dynamics function that describes the dynamics of its associated state variable over time. This modular representation enables independent analysis and optimization of individual components. In cases where R specifies purely mechanistic equations, the component dynamics are entirely defined using closed-form equations. Conversely, in a hybrid model, the mechanistic equation can be augmented with a neural network (implemented in PyTorch) to model residuals (i.e. in an 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 dynamics function as a separate code structure (for a generated HDTwin example, see Appendix I).Model optimization. The generated specification emits ω(θ), 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 ω * (θ) = arg min ω(θ) L(f θ,ω(θ) , D train ). In this work, we consider ω(θ) to be continuous variables, and as such, we optimize θ by stochastic gradient descent, using the Adam optimizer [32] . However, we note other optimization algorithms, such as black-box optimizers, could also be used (for more details, see Appendix F, Equation ( 5)). The parameter optimization step then yields the complete model, f θ,ω * (θ) .Quantitative evaluation. For each generated model, we evaluate them quantitatively. Specifically, we collect the validation mean squared error loss per component, which we denote as δ = [δ 1 , δ 2 , . . . , δ d X ] (Appendix F, Equation ( 6)). We also compute the validation loss of the overall model as well as υ = L(f θ,ω * (θ) , D val ). Finally, the generated model and its validation losses are included in a tuple and added to the top-K models P (g) ← P (g-1) ⊕ (f θ,ω * (θ) , δ, υ), where P (g) automatically removes the lowest performing models, and also only adds a new model to P (g) if it 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] .Model evaluation. The goal of the evaluation step is to reflect on the current set of top-K models, P (g) against requirements R and provide actionable and detailed feedback to the modeling agent for model improvement: H (g) ∼ LLM eval (R, P (g) ). We note that H (g) is provided in natural language and can be viewed as a 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 (1) engage in comparative analysis, identifying effective specifications in P (g) contributing to 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-leveraging the LLM's capacity to reason about proposed HDTwins to reflect these requirements via model improvement feedback.Enhancing search. By providing rich feedback to 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 componentspecific suggestions, proposing alternative decompositions, removing parameters, or introducing non-linear terms. It is noteworthy that the feedback H (g) , expressed flexibly in natural language, could easily be further enriched through direct human feedback. We demonstrate this human-inthe-loop capability by including expert feedback during the optimization process through H (g) 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 DTs. 
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 as a sequential modeling problem. In these settings, f θ,ω(θ) are typically black-box models, where θ ∈ Θ is the NN architecture and ω(θ) are its weights. Early 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 differential equations. These methods have found utility in modeling a 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 ([P1, 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 governing equations directly from data. Here θ ∈ Θ 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 ω(θ) from data. These techniques struggle to scale to higher-dimensional settings and rely on experts to 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 a 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. Correspondingly, 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 is an automated approach to jointly optimize hybrid model specification and its parameters. 
experiments and evaluation: 
 In this section, we evaluate HDTwinGen and verify that it significantly outperforms state-of-the-art methods in modeling system dynamics over time from an observed dataset and corresponding system description. 4Benchmark 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 human experts. Three are derived from a state-of-the-art biomedical Pharmacokinetic-Pharmacodynamic (PKPD) model of lung cancer tumor growth, used to simulate the combined effects of chemotherapy and radiotherapy in lung cancer [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 agent-based simulator (COVID-19) [65] , which is capable of modeling non-pharmaceutical interventions, such as physical distancing during a lockdown. Furthermore, we compare against an ecological model of a microcosm of algae, flagellate, and rotifer populations (Plankton Microcosm)-replicating an experimental three-species prey-predator system [66] . Moreover, we also compare against a 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 as D test , using the loss defined in Equation (5) and report this as T M SE . 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-art, we compare it with the most competitive and popular neural network models, which, when modeling the dynamics of a system over time, becomes a form of ODE model, that is a neural ODE [14] with action inputs (DyNODE) [68] . Moreover, we also compare against a recurrent neural network (RNN) [69] and a 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 a 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. 
main results: 
 We evaluated all our benchmark methods across all our datasets tabulated in Table 1 . HDTwinGen 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 This section provides an 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] out-of-distribution generalization, [P2] sample-efficient learning, and [P3] evolvability.[P1] Can an HDTwin generalize to out-of-distribution shifts? To explore out-of-distribution shifts, we adapt the Lung Cancer (with Chemo. & Radio.) to produce a training dataset of states in a range that is outside those observed in the test set over all trajectories (Appendix H.1). We tabulate this in Table 2 . Empirically, we find that HDTwinGen is more robust to out-of-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.[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, we 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 error over subsequent generations. We observe that HDTwinGen can indeed understand, reason, and iteratively evolve the generated code representation of the HDTwin to incorporate a 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 Number of Generations in accuracy and robustness. This iterative evolution process demonstrates HDTwinGen's 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 it to 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 an unobserved intervention of a 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 LLM and (2) adapted in 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 a 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 (Appendix H.5). Additionally, HDTwinGen is 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 H.10). 
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)-a promising approach that combines mechanistic 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 HDTwins consistently outperform existing approaches across multiple criteria, exhibiting superior out-of-distribution generalization, enhanced sample efficiency, and improved modular evolvability.Limitations. While our results are promising, several important limitations remain. HDTwinGen's efficacy depends critically on human experts providing initial system specifications and on the underlying LLM's domain knowledge and model generation capabilities. Our current implementation focuses exclusively on continuous-time systems, which, although broadly applicable, represent only a subset of real-world systems. Future work could extend our approach through human-in-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 a comprehensive evaluation of evolved models for fairness, bias, and privacy concerns before deployment in sensitive applications. We provide an illustrative example of HDTwinGen working in practice in Figure 5 .You are tasked with creating a PyTorch module to simulate the impact of chemo and radiotherapy on lung cancer tumor volume and chemotherapy drug concentration.Task Overview:• Objective: Simulate changes in tumor volume and drug concentration due to chemotherapy and radiotherapy dosages.• Key Variables: Tumor volume (cm3 ), drug concentration ( mg m3 ), chemotherapy dosage ( mg m3 ), and radiotherapy dosage (Gy). You are tasked with creating a PyTorch module to simulate the impact of chemo and radiotherapy on lung cancer tumor volume and chemotherapy drug concentration.Task Overview:• Objective: Simulate changes in tumor volume and drug concentration due to chemotherapy and radiotherapy dosages.• Key Variables: Tumor volume (cm3 ), drug concentration ( mg m3 ), chemotherapy dosage ( mg m3 ), and radiotherapy dosage (Gy). You are tasked with creating a PyTorch module to simulate the impact of chemo and radiotherapy on lung cancer tumor volume and chemotherapy drug concentration.Task Overview:• Objective: Simulate changes in tumor volume and drug concentration due to chemotherapy and radiotherapy dosages.• Key Variables: Tumor volume (cm3 ), drug concentration ( mg m3 ), chemotherapy dosage ( mg m3 ), and radiotherapy dosage (Gy). HDTwinGen can generate and further evolve HDTwins for a particular system based on user-given system requirements and a dataset D = {D train , D test } of state-action trajectories. First, the system requirements-which include dataset statistics are incorporated into a prompt and fed into the modeling agent that returns the code for the HDTwin. This HDTwin is then trained on the training dataset D train , and a validation loss is computed with D val . 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 on 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-performing HDTwin (w.r.t. validation performance) is returned. Overall, this produces an HDTwin that fulfills [P1-P3].Z A V V T 2 g n 9 v h t C A L b 2 1 X w R B i w V q E = " > A A A B 7 X i c b V B N S 8 N A E J 3 4 W e t X 1 a O X x S J 4 K k m h 6 r H g x W M F + w F t K J v N p l 2 7 y Y b d S a G U / g c v H h T x 6 v / x 5 r 9 x 2 + a g r Q 8 G H u / N M D M v S K U w 6 L r f z s b m 1 v b O b m G v u H 9 w e H R c O j l t G Z V p x p t M S a U 7 A T V c i o Q 3 U a D k n V R z G g e S t 4 P R 3 d x v j 7 k 2 Q i W P O E m 5 H 9 N B I i L B K F q p 1 R u H C k 2 / V H Y r 7 g J k n X g 5 K U O O R r / 0 1 Q s V y 2 K e I J P U m K 7 n p u h P q U b B J J 8 V e 5 n h K W U j O u B d S x M a c + N P F 9 f O y K V V Q h I p b S t B s l B / T 0 x p b M w k D m x n T H F o V r 2 5 + J / X z T C 6 9 a c i S T P k C V s u i j J J U J H 5 6 y Q U m j O U E 0 s o 0 8 L e S t i Q a s r Q B l S 0 I X i r L 6 + T V r X i X V d q D 9 V y v Z b H U Y B z u I At I = " > A A A C I H i c b V D L S g N B E J z 1 G e M r 6 t H L Y h A 8 h L A r k n g M 6 M F j B P O A b A i 9 k 0 k y Z H Z 2 m e k N h m U / x Y u / 4 s W D I n r T r 3 H y O G i S g o a i q p v u L j 8 S X K P j f F t r 6 x u b W 9 u Z n e z u 3 v 7 B Y e 7 o u K 7 D W F F W o 6 E I V d M H z Q S X r I Y c B W t G i k H g C 9 b w h z c T v z F i S v N Q P u A 4 Y u 0 A + p L 3 O A U 0 U i d X 9 h I v A B x Q E M l t 2 k k 8 Z I + Y o A I u 0 7 S w w h q B S F M v 7 e T y T t G Z w l 4 m 7 p z k y R z V T u 7 L 6 4 Y 0 D p h E K k D r l u t E 2 E 5 A I a e C p V k v 1 i w C O o Q + a x k q I W C 6 n U w f T O 1 z o 3 T t X q h M S b S n 6 t + J B A K t x 4 F v O i c X 6 0 V v I q 7 y W j H 2 r t s J l 1 G M T N L Z o l 4 s b A z t S V p 2 l y t G U Y w N A a q 4 u d W m A 1 B A 0 W S a N S G 4 i y 8 v k / p l 0 S 0 V S / d X + U p h H k e G n J I z c k F c U i Y V c k e q p EU = " > A A A B 8 n i c b V D L S g M x F L 1 T X 7 W + q i 7 d B I v Q V Z k R q S 4 L b l x W s Q + Y D i W T Z t r Q T D I k G a E M / Q w 3 L h R x 6 9 eQ Y z H g k W M Y G M l f x B j M y G Y Z w / z Y b X m N t w F 0 D r x C l K D A u 1 h 9 W s w k i S N q T C E Y 6 1 9 z 0 1 M k G F l G O F 0 X h m k m i a Y T P G Y + p Y K H F M d Z I v I c 3 R h l R G K p L J P G L R Q f 2 9 k O N Z 6 F o d 2 M o + o V 7 1 c / M / z U x P d B B k T S W q o I M u P o p Q j I 1 F + P x o x R Y n h M 0 s w U c x m R W S C F S b G t l S x J X i r J 6 + T 7 m X DSequence 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 model-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 a evolutionary multi-agent framework, while introducing more generalized mechanisms to incorporate loose-form prior knowledge.Discovering closed-form models. Closely aligned with our research are techniques aimed at discovering 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-dimensional settings and rely on 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 of prior knowledge at various stages of the search process [73, 74] . Furthermore, using LLMs to generate code, prior work has shown LLM multi-agent frameworks' ability to excel at large code-generation tasks [75, 76] , which we 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] .In the following, we outline the six real-world system dynamics datasets, where each dataset is either a real-world dataset or has been sampled from an accurate simulator designed by human experts.Three of our environments that we sample a dataset from are derived from a state-of-the-art biomedical Pharmacokinetic-Pharmacodynamic (PKPD) model of lung cancer tumor growth, used to simulate the combined effects of chemotherapy and radiotherapy in lung cancer [61] -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.)); for each model we sample a respective dataset. First, let us 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-the-art biomedical Pharmacokinetic-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 x(t) for days t after the cancer diagnosis-where the outcome is one-dimensional. The model has two binary treatments: (1) radiotherapy u r t and (2) chemotherapy u c t .EQUATIONWhere the parameters K, ρ, β c , α r , β r for each simulated patient are detailed in [61] , which are also described in Table 3 . Additionally, the chemotherapy drug concentration c(t) follows an exponential EQUATIONwhere the chemotherapy binary action represents increasing the c(t) concentration by 5.0mg/m 3 of Vinblastine given at time t. Whereas the radiotherapy concentration d(t) represents 2.0Gy fractions of radiotherapy given at timestep t, 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, p c and p r , are determined by the tumor diameter as follows:EQUATIONwhere D max = 13cm represents the largest tumor diameter, θ c = θ r = D max /2 and D(t) signifies the mean tumor diameter. The parameters γ c and γ r manage the extent of time-varying confounding.We useγ c = γ r = 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 x(0) ∼ U(0, 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 D train , D val , D test . Specifically for each benchmark method run for random seed, we re-sample the datasets. For each variation described above, we 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 as given in Appendix E.4.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 a population, and they can be in one of the following states minimally, of being susceptible to COVID-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 a simulation for a population of individuals. To ensure an accurate simulation, we simulate 24 countries collecting trajectories for each, wherein each simulation we use a population size of 1, 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, I(0) = U(10, 000, 100, 000), and forward simulate the simulation for 60 days. We repeat this process with independent random seeds to generate D train , D val , D test . 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.This describes an ecological model of a microcosm of algae, flagellate, and rotifer populations, thus replicating an experimental three-species prey-predator system [66] . We use the dataset made available by [79] 6 . The dataset consists of a single trajectory of 102 time steps, and we use a 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.This describes a 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 a 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.To be competitive we compare against popular black-box models, which, when modeling the dynamics of a system over time, becomes a form of ODE model, that is a 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 is a black-box neural network-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 a DyNODE model with 3-layer Multilayer Perceptron (MLP), with a hidden dimension of 128 units, with tanh activation functions, and make it competitive by using Xavier weight initialization [80] . 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 a 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.Causal Transformer is a state-of-the-art 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, respectivelywhich is unique to estimating counterfactual outcomes in treatment effect settings; we implemented only a single transformer to model the past outcomes, which is applicable to our datasets and task domains. Specifically, this consists of a 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 a linear layer, followed by the addition of a 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 a 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 a step learning rate scheduler of step size 1.0 and gamma 0.95; we also implement gradient clipping to 0.7, 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.Recurrent Neural Network [81] is a standard baseline that is widely used in autoregressive time series next step ahead prediction. We implement this where the input dataset is normalized to 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 to a linear layer to 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 a 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.Sparse Identification of Nonlinear Dynamics (SINDy) [11] , is a data-driven framework that aims to discover the governing dynamical system equations directly from time-series data, discovering a white-box closed-form mathematical model. The algorithm works by iteratively performing sparse regression on a library of candidate functions to identify the sparsest yet most accurate representation of the dynamical system.In our implementation, we use a polynomial library of order two, which is a feature library of L = {1, x 0 , x 1 , x 0 x 1 }. Finite difference approximations are used to compute time derivatives from the input time-series data, of order one. Here the alpha parameter is kept constant at 0.5 across all experiments, and the sparsity threshold is set to 0.02 for all experiments, apart from the COVID-19 dataset where it is set to 1 × 10 -5 .APHYNITY [60] is implemented using domain-specific expert models as defined in Appendix H.8 combined with a 3-layer MLP, with the same hyper-parameters as in [60] .Genetic programming (GP) is implemented using the implementation and hyper-parameters from the baseline in [82] .See the section Appendix E for the implementation details. Specifically, ZeroShot and ZeroOptim are 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-again using the same training as detailed in Appendix E.2.Our proposed method follows the framework as described in Section 4. We present pseudocode in Appendix E. H ∼ LLM eval (R, P (g) ) {Generate self-reflection on how to improve the HDTwin.} 10: end for 11: Return: f θ,ω(θ) * {The best fitting model that scored the lowest validation loss}Once the modeling agent has generated an HDTwin f θ,ω(θ) , 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 S context , of which examples of such a skeleton are given in Appendix E.4. However we stipulate that the skeleton must be a "torch.nn.Module", be called "StateDifferential", and the parameters must be initialized, and it must define a 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 LLM 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 a 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, to ensure fair comparison.Once the model is trained, we compute the val MSE and val MSE per component, which corresponds to the val loss per state output dimension Equation (6) . Notably, when we append the trained HDTwin back into P g 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.In the following we detail the prompt templates used. We always use the system prompt when interacting with the LLM. We consider the optimization loss of mean squared error on a dataset D and also consider a higher fidelity mean squared error loss per component. MSE Loss. Specifically, we optimize the following mean squared error objective,EQUATIONwhere N × T is the total number of state-action pairs in the dataset. For a given model find the parameters θ * that minimize this loss, i.e. θ * = arg min θ L(θ, D train ). Here we optimize θ by stochastic gradient descent, using the Adam optimizer [32] , however, we note other optimization algorithms could also be used such as 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 (j) to indicate the predictions for the j th component.ω m (θ * , D val ) = 1 N val × T Nval n=1 T n val i=0 f (j) θ,ω(θ) (x (n) (t i ), u (n) (t i ), t i )∆t -y (n) j (t i ) 2 (6)and collect these scalar validation losses per component into a vector ω = [ω 1 , ω 2 , . . . , ω m ], and compute its mean as v = 1 m m j=1 ω j (θ * , D val ) i.e. the validation loss.We employ mean squared error (MSE) to evaluate the benchmark methods on a held-out test dataset of state-action trajectories, denoted as D test , using the loss defined in Equation ( 5) and report this as T M SE . Each metric is averaged over ten runs with different random seeds, and we present these averages along with their 95% mean confidence intervals 7 . For each random seed run, we generate a new train, validation and test dataset independently, when we have access to a simulator. Additionally, when sampling a dataset from a 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 i9-12900K CPU @ 3.20GHz, 64GB RAM with an Nvidia RTX3090 GPU 24GB.To explore the out-of-distribution shifts, we adapt the Lung Cancer (with Chemo. & Radio.) simulator to have a training dataset of tumor volumes in a range that is outside of the tumor volume range within the test set over all trajectories. To do this we sampled a training dataset, with starting state tumor volumes sampled from x ∼ U(0, 574) and then collected 1, 000 trajectories for 60 time steps. We then sampled a test dataset out of distribution from that seen in training, with a starting initial state of x ∼ U(804, 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 ∆t = 1 24 rather than ∆t = 1 day resolution. We then also verified that the range of cancer volumes seen within the test set do not overlap at all with those in the training set, especially throughout and at the end of the trajectory.To assess whether HDTwinGen can understand and modify its HDTwin, we setup an experiment to see if the modeling agent can adapt an 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 a realistic intervention of a 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 β 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 a 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 is an SEIR model with black-box residual components, as seen in 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 as 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 learnt to optimize to.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 on average the HDTwin improves in the iteration loop in the beginning stages. 7 We use the code at https://stackoverflow.com/questions/15033511/ compute-a-confidence-interval-from-sample-data to compute these Here we list out the specific HDTwins that have been generated at generation steps and are labeled in Figure 3 . I 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 . We performed a complete re-run of our main experiments under the same settings, now using a 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 LLM model's performance as expected.Table 5 : Ablation of using different LLMs. Test MSE (T M SE ) averaged over ten random seeds.HDTwinGen is capable of using other LLM models, however, the best performance results are provided with better-performing LLMs (e.g. GPT-4). The results are presented with ± indicating 95% confidence intervals. We also explored the effect of changing the LLMs underlying temperature hyperparameter, here using GPT-4, varying the temperature from 0, 0.7 to 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 a different underlying temperature of the LLM.We conducted a complete re-run of our main experiments, ablating the prompt, which provides the prior information in the form of a textual prior to HDTwinGen. We provide the ablation results in Table 7 . Specifically the prompt as outlined in Appendix E is structured to include separately a {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 is also conveyed through the {skeleton code}, as this includes task-specific feature names as 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 a 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 taskspecific 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 HDTwinGen can still operate without any task-specific prior information, however having textual prior aids in generating better-performing models, and partially removing HDTwinGen operation instructions, makes it generate slightly less good models. 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 8 . Specifically, for COVID-19 modeling, we fit a SEIR model [84] , a Lotka-Volterra model for predator-prey 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.Modeling. The SEIR model is a compartmental model used in epidemiology to simulate how a disease spreads through a population. It 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: where N is the total population (assumed constant), β is the transmission rate, σ is the rate at which exposed individuals become infectious, and γ is the recovery rate. These parameters are crucial for capturing the dynamics of the disease spread and are estimated from data or literature.dS dt = -β SI N , dE dt = β SI N -σE, dI dt = σE -γI, dR dt = γI,Lotka-Volterra Model for Predator-Prey Dynamics. The Lotka-Volterra model describes the dynamics of biological systems in which two species interact, one as a predator and the other as prey.The model is represented by a set of first-order, non-linear, differential equations:dx dt = αx -βxy, dy dt = δxy -γy,where x and y represent the prey and predator populations, respectively. The parameters α, β, γ, 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-Volterra Triple Species Model. The extended Lotka-Volterra model incorporating a third species involves additional interactions that can represent various ecological relationships such as competition, predation, or mutualism. For the sake of illustration, let's consider a system with two predators and one prey. The model is described by the following set of differential equations:dx dt = x(α -βy -δz), dy dt = y(-γ + ϵx), dz dt = z(-µ + νx),where: x represents the prey population. y and z represent the two predator populations. α is the natural growth rate of the prey in the absence of predation. β and δ are the predation rates of the first and second predators on the prey, respectively. γ 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 = rN 1 - N K -C(N ) -R(N ),where N is the tumor cell population, r 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.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 a 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 of the underlying equation.Synthetic 1 (inc. γ sin(ωt)). Here the underlying equation isdx(t) dt = ρ log K x(t) -β c C(t) -(α r d(t) + β r d(t) 2 ) + γ sin(ωt) x(t) Synthetic 2 (inc. -δI(t)) dx(t) dt = ρ log K x(t) -β c C(t) -(α r d(t) + β r d(t) 2 ) -δI(t) x(t) Synthetic 3 (inc. log( K x(t)+N (t) )) dx(t) dt = ρ log K x(t) + N (t) -β c C(t) -(α r d(t) + β r d(t) 2 ) x(t) Synthetic 4 (inc. ϵ cos(ϕt)) dx(t) dt = ρ log K x(t) -β c C(t) -(α r d(t) + β r d(t) 2 ) + ϵ cos(ϕt) x(t) Synthetic 5 (inc. θC(t)d(t)) dx(t) dt = ρ log K x(t) -β c C(t) -(α r d(t) + β r d(t) 2 ) -θC(t)d(t) x(t)We observe in the tabulated results in Table 9 that HDTwinGen can still generate models that perform well. To investigate questions, of how well do the white-box models that HDTwinGen generates perform, we explore an ablation of HDTwinGen 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-Box)) it still performs acceptably, indicating that the white-box generated models are modeling well the underlying system when fitted to the dataset. Experts play an active role in model development in two main ways:1. Initial prompt. Experts can describe the system, and specify modeling goals and taskspecific requirements through S context (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 a fully white-box model through S context , 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 is the case and provide insight that the underlying LLM is able to interpret the HDTwin code model and selectively change parts. We provide a figure to illustrate this result, as seen in Figure 7 . Create a model to predict changes in tumor volume and chemotherapy drug concentration in NSCLC patients, using differential equations that incorporate logistic tumor growth with a log term, chemotherapy, and radiotherapy effects, with parameters estimated from patient data.Validation loss: 4.486Validation loss: 0.948 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 a fixed model budget. To assess this we performed the additional experiments of:1. Human-driven model development: Hiring two experienced software engineers (Mechanical 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 two experiments on the Lung Cancer (with Chemo. & Radio.) dataset, and the results are tabulated in Analyzing the results provides the following insights:1. 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 These are the final discovered hybrid models generated from our method HDTwin, for each respective environment. Re-e v a l u a t e t h e c o m p e t i t i o n t e r m s i n t h e e q u a t i o n s t o e n s u r e t h e y a r e b i o l o g i c a l l y m e a n i n g f u l and c o n s t r a i n t h e c o e f f i c i e n t s t o be non -n e g a t i v e d u r i n g o p t i m i z a t i o n . • 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. Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations. • The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?Answer: [NA] Justification: This paper does not release any pretrained models or collated datasets that might pose a risk to misuse.Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters. • Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images. • We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort.12. Licenses for existing assets Question: 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?Answer: [Yes]Justification: The paper cites original owners of all code (for baselines) and datasets used in the paper.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 a URL. • The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided. • If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a 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 creators. 
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/model as part of their submissions via structured templates. This includes details about training, license, limitations, etc. • The paper should discuss whether and how consent was 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. 
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 to participants and screenshots, if applicable, as well as details about compensation (if any)?Answer: [NA] Justification: This work does not involve crowdsourcing or 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 contribution of the paper involves human subjects, then as much detail as possible should be included in the main paper. • According to the NeurIPS 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. 
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 crowdsourcing nor research with human subjects.• 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.This model relates population size N , growth rate r, and carrying capacity K, ∂N /∂t = rN (1 -N /K)[30].Our implementation is available at https://github.com/samholt/HDTwinGen. We also provide a wider lab code repository at https://github.com/vanderschaarlab/HDTwinGen.COVASIM is an opensource simulator, from which we access it here https://github.com/ InstituteforDiseaseModeling/covasim.The Plankton Microcosm and Hare-Lynx datasets are both open source and available from https:// github.com/WillemBonnaffe/NODEBNGM.. ** R e f i n e I n t e r a c t i o n Terms * * : The i n t e r a c t i o n t e r m s c u r r e n t l y u s e s i m p l e m u l t i p l i c a t i v e f a c t o r s b e t w e e n p o p u l a t i o n s . C o n s i d e r more complex i n t e r a c t i o n s t h a t may i n c l u d e n o n l i n e a r e f f e c t s o r f u n c t i o n a l r e s p o n s e s , s u c h a s a H o l l i n g Type I I f u n c t i o n a l r e s p o n s e , t o b e t t e r c a p t u r e t h e d y n a m i c s o f p r e d a t i o n .. ** A d j u s t N e g a t i v e Growth R a t e * * : The o p t i m i z e d p a r a m e t e r ' a l p h a ' f o r t h e p r e y g r o w t h r a t e h a s a n e g a t i v e v a l u e , which i s b i o l o g i c a l l y u n r e a l i s t i c a s i t would i m p l y t h e p r e y p o p u l a t i o n d e c r e a s e s i n t h e a b s e n c e o f p r e d a t o r s . E n s u r e t h a t t h e g r o w t h r a t e p a r a m e t e r s a r e c o n s t r a i n e d t o be p o s i t i v e d u r i n g o p t i m i z a t i o n . 
 
 Appendix:  


Decision: 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.

Reviews: review: paper_summary: paper_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.
strength_weakness: strength_weakness: 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.

A 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:

- 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.
- 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.
- 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.
- 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.
questions: 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: Limitations clearly discussed.
paper_summary: paper_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.
strength_weakness: strength_weakness: 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. 
2. 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. 
3. Following-up from the previous point, more discussion on the computational complexity of the algorithm in the main paper could be beneficial. 
4. 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
5. 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.
questions: 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? 

In the limitations below, I have two questions that can be answered.
limitations: 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: 
- 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? 
- 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?  

Both points and questions above are particular relevant in critical domains and scenarios such as the medical scenarios used for demonstration in the paper. 

I 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.
paper_summary: paper_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.
strength_weakness: strength_weakness: 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. 

In 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). 

The 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. 

The 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.

Most of the experimental data seem to be low in dimension (at each time frame). Please clarify.
questions: 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.

It’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. 

For 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).


Clarification on the “spatial” dimension of the experimental datasets would be appreciated.
limitations: limitations: The authors provided adequate discussion about the limitation and future work for the presented work.
paper_summary: paper_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. 

Firstly, 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. 

Performance is reported on synthetic data from different models. Proposed model outperform the reported baselines.
strength_weakness: strength_weakness: - 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. 


- No code is provided.

- The work relies on a closed-sourced model (GPT-4) whose performance also depends on the specific checkpoints used.
questions: 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.

- 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.

- I'm unsure whether all the baselines are synthetic data or only some of them. Could you please clarify this?

- 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.

- Line 188: doesn't $P^{(g)}$ contain $w(\theta)$? If not, could you elaborate on why?

- Why are two LLMs needed? Couldn't a single one do both steps? Have you tried this?

- Suggestion: It'd be interesting to see how the algebraic expression evolve over time during training.
limitations: 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. 

- 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.

- 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.

- Besides, Table 10, the explainability potential of the model remains under-investigated.

Comments: title: Thank you
Thank you
Thank you
Thank you
comment: Thank you for your valuable feedback. We are glad to have addressed your concerns and appreciate your insights, which have significantly enhanced the quality of our work.
Thank you for the response. The new experiments seem to be informative. I would keep my originally positive score.
We appreciate your insightful feedback, which has been crucial in enhancing the quality of our submission. We are delighted to have successfully addressed your concerns and grateful for your constructive input throughout this process.
Thank you for the clarifications and the work put on getting the new baselines metrics. I think your work makes for a nice contribution, I will update the rate accordingly.
We're happy to have resolved your concerns and are thankful for your input, which has played a key role in refining the quality of our work.
Thank you for your valuable feedback. We are glad to have addressed your concerns and appreciate your insights, which have significantly enhanced the quality of our work.
Thanks for a thorough rebuttal which has clarified my main concerns -- the expansion of literature review and the addition of comparison to existing hybrid models is highly appreciated. I will raise my rating.
The authors have carefully considered my comments and improved the paper accordingly. As a consequence, I'm happy to increase my evaluation.