Learning The Delay in Delay Differential Equations

Authors: robert  stephany
maria  oprea
gabriella  torres
mark  walth
arnaldo  rodriguez-gonzalez
william  clark
Venue: ICLR_cc_2024_Workshop_AI4DiffEqtnsInSci
Published: 2024-02-10 02:23:25

Abstract: The intersection of machine learning and dynamical systems has generated considerable interest recently. Neural Ordinary Differential Equations (NODEs) represent a rich overlap between these fields. In this paper, we develop a continuous-time neural network approach based on Delay Differential Equations (DDEs). Our model uses the adjoint sensitivity method to learn the model parameters and delay directly from data. Our approach builds upon recent developments in NODEs and extends earlier neural DDE models, which assume the delay is known a priori. We rigorously justify our adjoint method and use numerical experiments to demonstrate our algorithm's ability to learn delays and parameters from data. Since the delay is rarely known \emph{a. priori}, our approach advances system identification of DDEs from real-world data.

Main paper: 
 introduction and motivation: 
 Neural ordinary differential equations (NODEs) have proven to be an efficient framework for various problems in machine learning Chen et al. (2018) . NODEs assume the map from input to target value can be modeled using a learned ODE. They solve an ODE for the forward pass and use an adjoint approach for backpropagation to update the network parameters. NODEs are a time and memoryefficient model for various regression and classification tasks Dupont et al. (2019) . We develop a novel continuous-time machine learning approach based on delay differential equations (DDEs).DDEs can model a wide range of phenomena: in computer networks, delays arise due to the transmission time of information packets Yu et al. (2004) ; in gene expression dynamics, a delay occurs due to the time taken for the messenger RNA to copy genetic code and transport macromolecules from the nucleus to the cytoplasm Busenberg & Mahaffy (1985) ; in population dynamics, the time taken for a species to reach reproductive maturity introduces a delay Kuang (1993) .There has been recent progress in developing Neural DDEs (NDDEs), the DDE counterpart of NODEs Anumasa & PK (2021) , Zhu et al. (2021 ), Zhu et al. (2022 ), Zhu et al. (2023) . For example, Zhu et al. (2021) and Zhu et al. (2022) show that NDDEs can learn models from time series data that NODEs struggle with. Our work extends Anumasa & PK (2021) , Zhu et al. (2021) , Zhu et al. (2022), and Zhu et al. (2023) , which assume the delay is known a priori (they can not learn the delay). By contrast, we derive an adjoint equation that allows our model to learn the delay, parameters, and initial conditions from data. This difference makes our approach more applicable as a system identification tool since the exact value of the delay is often unknown in practice.In this paper, we propose a novel algorithm to learn the delay and parameters of an unknown DDE using measurements of its solution. Our approach uses the adjoint sensitivity method (see Fig. 1 ), which we review in section 3.1. We implement our algorithm using PyTorch Paszke et al. (2019) . We then establish the efficacy of our approach on demonstrative examples in section 4. 
problem statement: 
 A DDE with a single delay, τ > 0, can be written in the form:EQUATIONwheref : R n × R n × R d → R n . Here, x 0 : [-τ, 0] → R n is the DDE's initial conditions * . Let C 2 ([0, T ], R n ) denote the set of C 2 functions from [0, T ] to R n .Our goal is to learn the parameters θ and delay τ from measurements of x. We do this using a loss functionJ : C 2 ([0, T ], R n ) → R +-with running and terminal cost, ℓ, g : R n → R + , respectively -of the form:EQUATIONWe consider the following scenario: We have N data measurements, {x(t j )} N data -1 j=0 , of the solution to an unknown DDE of the form of equation 1. Here, {t j } N Data -1 j=0 are known and partition [0, T ] (i.e., 0 = t 0 < t 1 < • • • < t N Data-1 = T ). Finally, we assume the unknown DDE has a constant initial condition given by x(t 0 ); i.e., x(t) = x(t 0 ) for t ∈ [-τ, 0] † . The main task of this paper is to solve the following minimization problem:Problem 1 Find arg min τ >0, θ∈R d J(x(•)) subject to ẋ(t) = f x(t), x(t -τ ), θ 0 < t ≤ T x(t) = x 0 , -τ ≤ t ≤ 0. 
methodology: 
 Unfortunately, solving problem 1 exactly is difficult. Therefore, we use an iterative approach based on gradient descent. We accomplish this by employing the adjoint sensitivity method to obtain the gradient of the loss with respect to the model parameters, θ, and the delay, τ .Assumption 1 In this paper, we assume that the running cost ℓ, the vector field f , and the terminal cost g are of class C 1 . We also assume that the map (θ, τ, t) → x θ,τ (t) (where x θ,τ is the solution of equation 1 -when the delay is τ and the parameters are θ -at time t) is of class C 2 .x(t)EQUATIONFigure 1 : Schematic drawing of the adjoint method for computing gradients. The dashed red curve on the upper half of the figure represents the computed solution to the original delay differential equation ẋ(t) = f (x(t), x(t-τ ), θ). We use the constant initial condition x(t) = x 0 for t ∈ [-τ, 0]. The dashed blue curve on the bottom half represents the computed solution to the adjoint problem.To derive our result, we use an approach similar to the one for ODEs in Rihan (2021) and Ayed et al. (2019) , but tailored for DDEs. We begin by defining the Lagrangian,L(x, λ) = J(x) + T 0 ⟨λ(t), ẋ(t) -f (x(t), x(t -τ ), θ)⟩dt.Here, λ : [0, T ] → R n is a Lagrange multiplier. We introduce λ due to the following observation:Remark 1 Let x : [-τ, T ] → R n be a solution to the DDE, equation 1. Then, for all continuous λ, L(x, λ) = J(x). Hence,∇ θ J(x) = ∇ θ L(x, λ) and ∂J ∂τ (x) = ∂L ∂τ (x, λ).Since J does not depend on λ, we can chose a λ that makes computing ∇ θ L and ∂L ∂τ more convenient.Theorem 1 Let ẋ = f (x(t), x(tτ ), θ) be a DDE with constant initial condition x 0 (t) = x 0 , t ∈ [-τ, 0] and define the cost as in equation 2. Assume that f , g, ℓ, and the flow x(t) satisfy Assumption 1. Choose λ : [0, T ] → R n to satisfy the following equation:     λ(t) = ∇ x ℓ (x(t)) -[∂ x f (x(t), x(t -τ ), θ)] T λ(t) -1 t<T -τ (t) [∂ y f (x(t + τ ), x(t), θ)] T λ (t + τ ) t ∈ [0, T ] λ(T ) = -∇ x g (x(T )) (3)Here, 1 t<T -τ is the indicator function on the set (-∞, Tτ ). Likewise, ∂ x f , ∂ y f , and ∂ θ f denote the (Frechet) derivative of f with respect to its first, second, and third arguments, respectively. Then, the derivatives of the cost function with respect to θ, τ and x 0 are:EQUATIONWe prove this result in Appendix B.We train a model to learn the delay and parameters of the DDE that generated the data. Algorithm 1 summarizes our approach. We begin with random initializations of θ and τ . Each epoch, we solve Algorithm 1: Learning delay and parameters from dataInput: {x(t j )} N data -1 j=0 Output: θ, τ 1 Initialization Set x 0 = x(t 0 ). Select random, initial values for θ, τ . 2 for i ∈ {1, . . . , N epochs } do 3 Solve ẋ(t) = f (x(t), x(t -τ ), θ) using the initial condition x(t) = x 0 , t ∈ [-τ, 0]. 4Compute J(x) and check for convergence. Break if converged.Otherwise, solve for the adjoint, equation 3 6Compute ∇ θ J(x) and ∂J(x)/∂τ using equation 4.Update θ and τ using ADAM. 8 end for equation 1 and compute the loss using the model's current θ and τ values. We then use this trajectory to solve equation 3 backward in time to find the adjoint. Next, we use equation 4 to compute the gradient of the loss. Finally, we use the Adam optimizer to update θ and τ Kingma & Ba (2014) . We repeat this for N epochs epochs or until the loss drops below a user-specified threshold.We implement Algorithm 1 using PyTorch Paszke et al. (2019) . We wrap the forward and backward procedures in a torch.autograd.Function class. In the forward pass, we use a secondorder Runge-Kutta solver to calculate the predicted trajectory. In the backward pass, we use the same DDE solver to calculate the adjoint by solving equation 3 backward in time. Our code then computes and returns the gradient of the loss with respect to τ and θ.In our implementation, we approximate the integral in the loss, equation 2, using a trapezoidal rule. We use SciPy's interpolate function to interpolate between the data points, {x(t j )} N data -1 j=0 . We then evaluate this interpolation and the predicted trajectory at the quadrature points to compute the quadrature rule. Our implementation is open-source and is available at https://github. com/punkduckable/NDDE. 
experiments and numerical results: 
 Here, we test Algorithm 1 by learning the parameters and the delay in a Logistic Delay Equation model. We repeat this analysis on a Delay Exponential Decay model in Appendix A. For these experiments, we use T = 10, g(x) = ∥x(T )x(T )∥ 2 2 , and ℓ(x(t)) = ∥x(t)x(t)∥ 2 2 , where x denotes the true solution to the hidden DDE (the data set consists of measurements of x). Further, we generated all our datasets by solving the corresponding equation using the forward Euler method with a step size of dt = 0.1. We use the resulting discretized solution (with 100 data points) as the target dataset, {x(t j )} N Data -1 j=1 , t j = j * 0.1. Finally, we use the Adam optimizer with a learning rate of 0.1, β 1 = 0.9, and β 2 = 0.999 Kingma & Ba (2014) .The Logistic Delay Equation was first proposed to understand oscillatory phenomena in ecology. It can model the dynamics of a single population growing toward saturation Ruth & Gergely (2020) . For this experiment, the true DDE is:ẋ(t) = x(t) (1 -x(t -1)) , t ∈ [0, T ] x(t) = 2 t ∈ [-τ, 0]. We use the model f (x, y, (θ 1 , θ 2 )) = θ 0 x(1 -θ 1 y).We initialize our model with θ 0 = θ 1 = 1.75 and τ = 1.75. We then use Algorithm 1 to learn θ and τ . We use N epochs = 1, 000 and set the loss threshold at 0.001. The loss drops below this threshold after just 144 epochs. Table 1 reports the discovered values of θ 0 , θ 1 , and τ . Thus, our approach identifies the correct parameters and τ values. Fig. 2b and 2a depict the loss landscape for this model; they show that the loss has a sink around the true θ and τ values. Fig. 2a also exhibits the path our algorithm takes as it trains. Finally, Fig. 3b shows that the target and final predicted trajectories are closely aligned.  
conclusions: 
 In this paper, we presented an algorithm to learn the time delay, parameters, and initial condition from data. Our approach employs the adjoint sensitivity method, suitably adapted for DDEs.We provided a rigorous derivation of our methodology and validated the efficacy of our approach through numerical experiments. In practice, the delay is rarely known a priori. Therefore, our contribution represents a crucial advance in system identification for real-world systems governed by delay differential equations.There are many potential future research directions to extend our results. We would first like to modify our approach to work with arbitrary initial conditions rather than constant functions. Doing so requires introducing a second Lagrange multiplier, µ, which subsequently changes the adjoint equations, equation 3. Second, from a numerical perspective, we want to understand how robust our algorithm is to limited and noisy data. We would also like to understand how various initialization schemes for θ and τ impact this robustness. Lastly, we want to develop theoretical error bounds for the discovered delay and a generalization error for the learned model. Delay Exponential Decay Equation: The Delay Exponential Decay Equation arises in the study of cell growth, specifically when the cells need to reach maturity before reproducing Rihan & Bocharov (200) . For this experiment, the true DDE is:ẋ(t) = -2x(t) -2x(t -1) t ∈ [0, T ] x(t) = 2 t ∈ [-τ, 0]For this experiment, we use the model:f (x, y, θ) = θ 0 x + θ 1 y.We initialize our model with θ 0 = θ 1 = -3.0 and τ = 2.0. We then use Algorithm 1 to learn θ and τ . Here, we use N epochs = 1, 000 and a loss threshold of 0.001. The loss drops below this threshold after 505 epochs. Table 2 reports the discovered values of θ 0 , θ 1 , and τ . Fig. 3a shows the final predicted and target trajectories. Thus, our approach correctly identifies the model parameters and can faithfully reproduce the target trajectory. 
 
 Appendix: In this appendix, we use an adjoint approach to derive an expression for the gradient of J with respect to λ, τ , and x 0 . ‡ First, we will consider the gradient with respect to θ. Fix (θ, τ, x 0 ) ∈ R d × (0, ∞) × R n . Since ∇ θ J(x) = d i=1(∂J(x)/∂θ i ) e i , it suffices to consider the case when θ ∈ R; we can extend to θ ∈ R d case by applying R case component-wise. Throughout this proof, we use results on Frechet differentiation; see Munkres (2018) or an equivalent reference.Before we proceed, we need some assumptions and notation: To make clear that the solution x depends on θ, we write x θ in place of x. Further, we assume that f , g, and ℓ are of class C 1 (each component function has continuous partial derivatives). Finally, we assume the map (θ, t) → x θ (t) is of class C 2 (each component function has continuous second-order partial derivatives).Next, let us clarify our derivative notation. This proof involves derivatives of functions of the form a : R n → R m , b : R n → R, and c : R → R. We use the notation ∂ x a(x) ∈ R m×n to denote the matrix of the (Frechet) derivative of a at x. If n = 1, we treat the m × 1 matrix ∂ x a(x) as a vector in R m . For the second class, we use the notation ∇ x b(y) ∈ R n to denote the gradient of b at y ∈ R n . Finally, for the third class, we use the notation ∂c(z)/∂z ∈ R to denote the (partial) derivative of c with respect to z. Time derivatives are the one exception to this notation: we write ḣ(t) for the time derivative h, regardless of whether h is real, vector, or matrix-valued.With that cleared up, let us consider ∂J(x θ )/∂θ. Doing so reveals a problem, which subsequently motivates the Lagrangian. Since the maps θ → x θ (T ) and u → g(u) are C 1 , they are differentiable. Therefore, by the chain rule, the map θ → g (x θ (T )) is differentiable with gradient⟨∇ x g (x θ (T )) , ∂ θ x θ (T )⟩. Similarly, for each t ∈ [0, T ], the map θ → ℓ(x θ (t)) is differentiable with gradient ⟨∇ x ℓ(x θ (t)), ∂ θ x θ (t)⟩.Critically, the components of the latter gradient are continuous on [0, T ]. Therefore, we are justified in differentiating under the integral sign. In particular, the map θ → T 0 ℓ(x θ (t)) dt must be differentiable with gradientT 0 ⟨∇ x ℓ(x θ (t)), ∂ θ x θ (t)⟩ dt. Therefore, the map θ → J(x θ ) is differentiable with ∂J ∂θ (x θ ) = ⟨∇ x g (x θ (T )) , ∂ θ x θ (T )⟩ + T 0 ⟨∇ x ℓ (x θ (t)) , ∂ θ x θ (t)⟩ dt.Unfortunately, evaluating this expression requires knowing ∂ θ x θ (t) for each t ∈ [0, T ]. Deriving an expression for this quantity is possible but non-trivial. Thankfully, we can address this problem by adding zero to J in a clever way: this engenders the LagrangianEQUATIONHere, λ θ : [0, T ] → R n is a Lagrange multiplier. As we stated in remark 1, λ θ can be any continuous function. However, we restrict ourselves to λ θ for which the map (θ, t) → λ θ (t) is of class C 1 .We can express the Lagrangian as a single integral using the formula for J:EQUATIONWe prove that L is differentiable with respect to θ and derive a formula for the corresponding gradient. We already know that the maps θ → T 0 ℓ (x θ (t)) dt and θ → g (x θ (T )) are differentiable. Now, fix t ∈ [0, T ] and consider the map θ → ⟨λ θ (t), ẋθ (t)f (x θ (t), x θ (tτ ), θ)⟩. To begin, by assumption, the map θ → x θ (t) is of class C 2 . In particular, this means that the component functions of the map θ → (x θ (t), x θ (tτ ), θ) ∈ R n × R n × R have continuous partials. Thus, this map must be differentiable. Since we assume that f is also of class C 1 , it must be differentiable. Therefore, by the chain rule, the map θ→ f (x θ (t), x θ (t -τ ), θ) is differentiable with ∂ θ f (x θ (t), x θ (t -τ ), θ) = (∂ x f (t)) ∂ θ x θ (t) + (∂ y f (t)) ∂ θ x θ (t -τ ) + ∂ θ f (t).For brevity, let f (t) be an abbreviation of f (x θ (t), x θ (tτ ), θ). Since x θ is of class C 2 , the map θ → ẋθ (t) must be differentiable as well. Further, since (t, θ) → x θ (t) is C 2 , the equality of mixed partials tells us that∂ θ ẋθ (t) = ∂θ x θ (t).Finally, since λ θ is of class C 1 , the component functions of both arguments of the inner product ⟨λ θ (t), ẋθ (t)f (x θ (t), x θ (tτ ), θ)⟩ have continuous partial derivatives with respect to θ. Therefore, the map θ → ⟨λ θ (t), ẋθ (t)f (x θ (t), x θ (tτ ), θ)⟩ must be differentiable with∂ ∂θ λ θ (t), ẋθ (t)-f (t) = ∂ θ λ θ (t), ẋθ (t) -f (x θ (t), x θ (t -τ ), θ) + λ θ (t), ∂θ x θ (t) -(∂ x f (t)) ∂ θ x θ (t) -(∂ y f (t)) ∂ θ x θ (t -τ ) -∂ θ f (t) .Crucially, since ẋθ (t) = f (x θ (t), x θ (tτ ), θ), the first term in this expression vanishes. Thus,EQUATIONBy inspection, each component of each argument of the inner product on the right side of equation 7 is continuous. Therefore, we are justified in differentiating under the integral. In particular, the map θ → L(x θ , λ θ ) must be differentiable with gradient∂ ∂θ L(x θ , λ θ ) = T 0 ⟨∇ x ℓ(x θ (t)), ∂ θ x θ (t)⟩ + λ θ (t), ∂θ x θ (t) -⟨λ θ (t), (∂ x f (t)) ∂ θ x θ (t) + (∂ y f (t)) ∂ θ x θ (t -τ ) + ∂ θ f (t)⟩ dt (8) + ⟨∇ x g (x θ (T )) , ∂ θ x θ (T )⟩ .To proceed, we need to make a few simplifications to the expression above. We will focus on two parts of the integrand, which we colored teal and violet. We begin with the teal portion of the integrand. Using integration by parts,EQUATIONTo get the second equality, note that x θ (0) = x 0 . Since x 0 is a constant, ∂ θ x θ (0) = 0 ∈ R n . We can now move to the violet portion of the integrand. In particular,EQUATION)To get from the second to the third line, we observe that for t < τ , x θ (tτ ) = x 0 , which means that ∂ θ x θ (t) = 0. Thus, the blue integral is zero. In the red integral, we redefine t as tτ . In equation 10, ∂ ∂θ L(x θ , λ θ ) = T 0 ⟨∇ x ℓ(x θ (t)), ∂ θ x θ (t)⟩ -λθ (t), ∂ θ x θ (t) -⟨λ θ (t), (∂ x f (t)) ∂ θ x θ (t)⟩ + 1 t<T -τ (t) ⟨λ θ (t + τ ), (∂ y f (t + τ )) ∂ θ x θ (t)⟩ + ⟨λ θ (t), ∂ θ f (t)⟩ dt + ⟨∇ x g (x θ (T )) + λ θ (T ), ∂ θ x θ (T )⟩ = T 0 ∇ x ℓ (x θ (t)) -λθ (t) -[∂ x f (t)] T λ θ (t)- 1 t<T -τ (t) [f y (t + τ )] T λ θ (t + τ ), ∂ θ x θ (t) dt (11) + ∇ x g (x θ (T )) + λ θ (T ), ∂ θ x θ (T ) + T 0 ⟨λ θ (t), ∂ θ f (t)⟩ dt.Thus far, our only assumption about λ θ is that the map (t, θ) → λ θ (t) is C 1 . Hence, we are free to choose any λ θ . In particular, we can select λ θ to satisfy the following differential equation:EQUATIONWe refer to this choice of λ θ as the adjoint. Substituting the adjoint into equation 11 givesEQUATIONTo generalize this to a vector valued θ, notice that if we replace θ with θ i in equation 13, the integral on the right is the i'th component ofEQUATIONWith this established, let us now turn our attention to τ . As we did with θ, we will fix (θ, τ, x 0 ) ∈ R d × (0, ∞) × R n . We will also denote the solution, x, by x τ to make its dependence on τ explicit.Using the same argument we used in the θ case, we introduce a Lagrange multiplier λ τ : [0, T ] → R n with the restriction that the map (τ, t) → λ τ (t) is of class C 1 :EQUATIONWe now argue that L is differentiable with respect to τ . We can reuse most of our argument from the θ case with one caveat: we need to re-established the differentiability of the map τ → f (x τ (t), x τ (tτ ), θ) (for a fixed t). Before, we proved this by showing that each of f 's arguments is differentiable. By assumption, the maps τ → x τ (t) and τ → θ (the latter being a constant) are differentiable. However, the second argument, x τ (tτ ), depends on τ both explicitly (in its argument, tτ ) and implicitly (because changing τ changes the DDE solution). However, this map is the composition of the map τ → (τ, t-τ ) with the map that sends (τ, t) to the solution of the DDE at time t when the delay is τ . By assumption, both maps are C 1 and thus differentiable. Hence, by the chain rule, the map τ → x τ (tτ ) is differentiable with derivative ∂ τ x τ (tτ )ẋτ (tτ ).Having established this map is differentiable, we can conclude (again by the chain rule) that the map τ → f (x τ (t), x τ (tτ ), θ) is differentiable. We can now proceed exactly as in the θ case to conclude that L is a differentiable function of τ . In particular, We can treat the teal and violet portions of the integrand as we did before. Further, since x τ (t-τ ) = x 0 for t < τ , the orange portion of the integrand is zero in [0, τ ]. Therefore,EQUATION+ ∇ x g (x τ (T )) + λ τ (T ), ∂ τ x τ (T ) + We can use a nearly identical argument to prove that J is a differentiable function of x 0 . Specifically,EQUATIONHere, λ x0 also satisfies equation 12. Dropping the θ, τ , and x 0 subscripts in equation 14, equation 18, and equation 19, respectively, gives equation 4. □‡ In our experiments, for simplicity, we assume that x0 is known. However, this assumption is not necessary as it is possible to learn this parameter from data. See equation 4. 


Decision: Dear Authors,

Thank you for submitting the draft.

Both reviewers agree that the presented work presents interesting strengths. It is expected that authors will be addressing comments by the reviewers in the final draft. 

regards

AC

Reviews: review: title: title: Review of "Learning The Delay in Delay Differential Equations"
main_review: main_review: # Summary

The paper derives an adjoint method that allows the delay in neural delay differential equations to be learned rather than specified *a priori*.
In experiments on simple systems, the proposed algorithm is able to recover the correct parameters and delay from the training data.

# Strengths
1. (**significance**)
The paper builds well upon previous work and addresses an interesting, impactful problem.
2. (**quality**)
The main contribution of the paper is the proof of the adjoint equation in Appendix B, which is thorough and well written, although I haven't verified its correctness.
3. (**clarity**)
More generally, the presentation of the paper, including figures, theorems, and algorithms, is excellent.
4. (**completeness**)
The authors set out extensive and clear opportunities for further work.

# Weaknesses
1. (**completeness**) 
While the question of "limited and noisy data" is mentioned briefly in the conclusion, I believe this is a fairly fundamental issue for the method that deserves detailed discussion.
In particular, the method requires interpolating between data points, since the delay $\tau$ is treated as a continuous variable while the data is necessarily discrete.
I strongly suspect that the interpolation may struggle in the presence of noisy, sparse, or irregularly sampled data.
Therefore, a complete discussion of the limitations of this method should include a detailed examination of this point.
1. (**completeness**) 
The experiments, while a reasonable proof of concept, are perhaps somewhat trivial, consisting each of only two learnable parameters in addition to the delay.
I would be very interested to know how well the proposed algorithm performs in the presence of many more parameters; for example, in [1] the authors train NDDEs with $\mathcal{O}(100k)$ parameters.
1. (**completeness**) 
Another discussion missing from the paper is the possibility of using a *discretize-then-optimize* approach to training NDDEs, i.e., applying automatic differentiation directly to the operations of the numerical ODE solver, as opposed to *optimize-then-discretize*, i.e., adjoint sensitivity analysis; see [2].
While I'm not certain of this, it may be possible to skip step 5 of Algorithm 1 entirely and do step 6 using automatic differentiation instead of using Eq. 4.
At the very least, I think it would be interesting to compare such an approach to the adjoint method derived in this paper.
1. (**completeness**) 
There appears to be some closely related work that is not mentioned in the paper [3].

# Other Comments
1. A link to a GitHub repository is included on page 3, but I didn't click on it and I don't think any identifying information about the authors can be determined from the URL.
1. One easy-to-fix problem with the paper's presentation is the bibliography, where some references are given to arXiv preprints instead of the relevant conference proceedings.

# Conclusion
Overall, I like this paper.
The weaknesses discussed above are intended only as constructive criticism that might benefit future versions of this work.

# Citations
[1] Qunxi Zhu, Yao Guo, and Wei Lin. Neural Delay Differential Equations. In *International Conference on Learning Representations*, 2021.

[2] Patrick Kidger. On neural differential equations. arXiv preprint arXiv:2202.02435, 2022.

[3] Xunbi A. Ji, Gábor Orosz. Learning Time Delay Systems with Neural Ordinary Differential Equations. IFAC-PapersOnLine, Volume 55, Issue 36, 2022.
title: title: Important contribution to extend neural ODEs to delay differential equations
main_review: main_review: The authors show a new method learn the delay parameter alongside 'normal' parameters of neural delay differential equations by proposing  a new algorithm to calculate the adjoint for neural delay differential equations. With two simple examples, the authors prove that their algorithm can retrieve the data generating delay parameter and normal parameters.

This work could be a great starting point for developing neural delay differential equations further. The authors give some examples for systems governed by delay differential equations and I agree that the topic is highly relevant. Here, the relevance could have highlighted better by finding prominent examples from different domains. The work is original and can advance the use of neural differential equations to new areas by allowing for delay parameters. A point that may be missing is an ablation experiment that checks if a delay of 0 is found when the neural delay differential equations is confronted with data generated with an ODE. 

Pro:
- extends neural ODEs to delay differential equations
- seems to be the first time that adjoint for neural delay differential equations is described
- delay parameter can for the first time be estimated in neural delay differential equations

Con:
- some sloppiness with respect to formatting citations - (Author, Year) formatting is not used when appropriate. This makes reading a bit tedious at times.