HyperSINDy: Deep Generative Modeling of Nonlinear Stochastic Governing Equations

Weaknesses

1. The problem setting is confusing and the baseline models are not adequate. Details see the questions part.

Response: Please see responses below.

Questions

1. Problem setting about the stochastic system. As the author mentioned, they used an alternative definition of the stochastic dynamics, e.g random ODE. However, the reviewer is confused by the benefit of such definition compared to the deterministic setting. In result part 4.1, both the mean and the STD of the underlying system are shown. However, the STD form doesn't corresponds to any dynamics and the mean is close the the true mean of the system. In such circumstances, it seems only the mean estimation is important and the std cannot be leveraged to judge the performance of the proposed model. For part 4.2, the STD result is also confusing. It is compared with the diffusion terms but it is totally different from the true diffusion term. Therefore, the reviewer wants to ask why we need to include the STD and how it can help the proposed model.

Response: We apologize for this confusion, which we believe stems from the interpretation of the RDE used to simulate the dynamics. In each case, mean and STD refer to the parameters of the (Gaussian) distribution associated with each coefficient; trajectories are simulated by constructing sample paths from these distributions. In other words, each coefficient is effectively coupled to an independent Wiener process; the system is parameterized by an n-dimensional Wiener process. Accordingly, each ground truth trajectory is associated with only one particular sample path for the n-dimensional Wiener process.

The STDs factor into the dynamics by specifying the level of dynamical noise that drives each coefficient. The relevance of the STD to the dynamics is most clear in the rows of Fig. 2 corresponding to $\sigma=10$. In these cases, the ground truth trajectory (“Original”) is highly stochastic. HyperSINDy succeeds in recovering the ground truth equation along with the correct coefficient means. However, the mean (middle column) is a poor representation of the original system, which behaves much more randomly. Thus, by discovering appropriately large STD values, HyperSINDy is able to accurately represent the level of uncertainty in the system’s evolution; this also enables HyperSINDy to faithfully recapitulate the dynamical behavior of the original system. E-SINDy fails to capture the appropriate level of uncertainty (as it assumes no dynamical noise); as such, it cannot faithfully recapitulate the dynamical behavior of the original system (Fig. S4).

2. The result needs to compare with more SOTA models and include more comprehensive metrics. There are several model combining learning based method and SINDy-like algorithm for the equation discovery tasks [1][2]. Also, by checking the Figure 2, we could find that the discovery form is different under different sigma . For a equation discovery problem, it is important to get a consistent and correct form. Therefore the reviewer suggests adding more metrics on evaluating if the proposed model can get the correct form, e.g, precision and recall metrics.

Response: We believe that this point is partially addressed by our clarification of the current problem setting as stochastic rather than merely noisy, deterministic dynamics. Nonetheless, we completely agree that comparison against another state-of-the-art method for probabilistic equation discovery, such as Bayesian Spline Learning, would improve the manuscript. We have updated the manuscript to include comparisons with this method along with additional quantitative metrics – thank you for these suggestions. Refer to Table S4 in the appendix for precision and recall results that complement the experiments ran for Table 1 in the main text (on Lorenz and Rossler systems generated with 10 seeds). Refer to Table S5 in the appendix for RMSE, precision, and recall comparisons between HyperSINDy and E-SINDy on multiple 10D Lorenz-96 trajectories, each generated with a different noise level. Refer to Table S6 in the appendix for RMSE, precision, and recall comparisons between HyperSINDy, E-SINDy, and Bayesian Spline Learning on multiple Lotka-Voltera trajectories simulated as RDEs with Gaussian noise on every coefficient. (response continued below)

I have changed my score to reflect these efforts.

Weakness

1. The capacity of (theta(x) still holds as a constraint for the performance, especially in the high-dimensional setup. It would be great if the authors could discuss the impact of theta(x) . For example, what would the performance be if certain symbolic terms (shown up in the true equations) are missing in the dictionary theta(x).

Response: We agree that the model is constrained by the capacity of the function library. This is true of most equation discovery methods. The polynomial basis has broad applicability; moreover, there is often a priori knowledge motivating choice of specific classes of basis functions. The SINDy framework has had widespread success with alternative libraries (e.g., Kaptanoglu et al., 2022). Our method is amenable to such alternative libraries, and thus retains the flexibility of the SINDy framework.

Questions

1. What is the column of ''STD'' in Figure 2 showing? Are they showing the standard deviation of the estimates? If that is the case, plugging in the standard deviation as the coefficients in the equations are confusing.

Response: Thank you for pointing out this confusion. Correct, this figure is showing the standard deviation of the estimated coefficients, although these standard deviations correspond to the noise processes that contribute to the dynamics. We have updated the figure to clarify the interpretation of the standard deviations.

2. It would be great if the authors could provide more evaluation metrics for generated trajectories. Metrics like Lyapunov exponents would be helpful to see how good the performance is.

Response: Thank you for this excellent suggestion. Although we were unable to carry out a systematic analysis of the Lyupanov exponents for generated trajectories, we have incorporated two additional metrics, precision and recall (Tables S4-S6). These additional metrics enable assessment of the ability to recover the correct terms of the original equation.

3. How robust the performance would be across different choice of the dimension of z?

Response: Thank you for raising this important question. We have carried out additional experiments to assess performance as a function of latent dimension. Specifically, we train HyperSINDy models with varying dimensions of z on each of the 3D Stochastic Lorenz trajectories ($\sigma = 1$, $\sigma = 5$, $\sigma = 10$).

Results from this procedure are shown in Fig. S6. We plot the RMSE of the mean and standard deviation of discovered coefficients, as compared to ground truth. It is true RMSE can vary as a function of z dimension. It is noteworthy, however, that for each trajectory, RMSE is always lower than E-SINDy, regardless of z dimension (refer to Table 1 for E-SINDy RMSE). These results help validate our choice of z in the paper: we use 2x the state dimension for z, and this seems to have good performance at balancing RMSE of both the mean and standard deviation of discovered coefficients.

References

Kaptanoglu, A.A. et al. (2022) ‘PySINDy: A comprehensive Python package for robust sparse system identification’, Journal of Open Source Software, 7(69), p. 3994. Available at: https://doi.org/10.21105/joss.03994.

Thank you for your prompt follow-up. As elaborated in our latest response to Reviewer 7UM6, we stress that our manuscript does not claim the ability to explicitly map between particular RDE and SDE expressions, and our results are not contingent on the ability to carry out this transformation. Rather, we simply note that the conjugacy between the two formulations means that the results obtained in one framework are recognized to be directly pertinent to those obtained in the other (Han and Kloeden, 2017). Thus, our work retains relevance to a very broad class of problems that have been modeled as SDEs – ultimately, SDEs and RDEs are simply two alternative modeling choices for studying dynamics whose evolution involves some uncertainty.

Note that the Lotka-Volterra SDE problem (Fig. 3) demonstrates the utility of our approach for a problem explicitly formulated as an SDE -- HyperSINDy captures the appropriate dynamics (as confirmed by evaluating the Kramers-Moyal coefficients) and, importantly, identifies the physical terms that are present in the drift component of the original SDE, thus supporting insights into the physical interactions in the system.

We will update the Appendix to include these clarifications, as well as an example transformation for multiplicative noises (see latest response to Reviewer 7UM6.

Weaknesses

1. Experiments are conducted in simulated environments where the simulation parameters match to the modeling assumptions (mainly around Gaussianity). I would love to see more experiments confirming the applicability of the approach to broader problems, especially with real data.

Response: The insight exploited in this paper is that the standard VAE gaussian prior is well-matched to the well-studied case of stochastic dynamics involving white noise. We completely agree that extension to other kinds of noise will be an important direction for future extensions. As an initial investigation into the feasibility for such extensions, we have included a new experiment (Fig. S5) that supports the potential for HyperSINDy to model more complicated noise types (here, noise drawn from a half-Normal distribution).

Questions

1. As mentioned above, all the experiments are conducted using Gaussian distributions which match to the posterior distribution assumed for variational inference. Can authors comment on the limitations of these experiments?

Response: (see above response)

2. The approach aims to learn both the functional form and parameters of the differential equations. Even though I agree that this might help with interpretability, I worry that the identifiability issues might be prominent. Do the authors expect any identifiability problems?

3. The promise of learning functional form is achieved through the function library. Are there any limitations of using such an approach?

Response: We agree that there can be identifiability issues in equation learning; addressing these issues generally comes at the cost of reduced flexibility. This is precisely the tradeoff reflected in these questions. Thus, specifying a function library and imposing sparsity are significant aids toward identifying a parsimonious model. Moreover, by assuming the structure of the noise, we avoid the major identifiability issues that arise when attempting to learn the form of both the deterministic and stochastic terms. These constraints each necessitate a loss of generality, though they are appealing as they retain very broad applicability. The primary limitation of a function library is that it requires a priori knowledge about the dynamics, i.e. some inductive bias is necessary. However, this inductive bias is also the primary motivation behind equation discovery methods: we want to discover a sparse equation from a finite set of possible terms that could describe the system’s dynamics. One alternative might be to enable a more fully data-driven search for the basis functions (e.g., use via a genetic algorithm (Schmidt and Lipson, 2009)), although such approaches are generally incur much greater computational cost.

4. What are the limitations of using a Gaussian prior with diagonal covariance for the generative model?

Response: One limitation is that such a prior assumes the true latent vector is built of multiple i.i.d. random variables (diagonal covariance). However, one could certainly imagine noise processes for which this assumption does not hold. Even though the hypernetwork contains nonlinearities (and thus can model dependencies between coefficients), some systems may still be challenging to learn. Future work could explore using alternate priors, such as Gaussian Process Variational Autoencoder (Casale et al., 2018), which may better account for such dependencies.

Despite this limitation, we note that the network is able to learn more complex distributions that include dependencies between coefficients. In particular, refer to Figure S5 in the appendix for results on a Lotka-Volterra system simulated with coefficients drawn from a Half-Normal distribution. Even though HyperSINDy uses a standard gaussian prior with diagonal covariance, it is able to approximate the shape of the distribution of the true non-Gaussian coefficients much better than E-SINDy. We view this experiment as a proof of principle for the flexibility provided by hypernetwork’s nonlinearities: it is able to learn a mapping from latent space (a simple standard gaussian) to a much more complex coefficient space.

References

Schmidt, M. and Lipson, H. (2009) ‘Distilling Free-Form Natural Laws from Experimental Data’, Science, 324(5923), pp. 81–85. Available at: https://doi.org/10.1126/science.1165893.

Casale, F.P. et al. (2018) ‘Gaussian Process Prior Variational Autoencoders’, in Advances in Neural Information Processing Systems. Curran Associates, Inc. Available at: https://proceedings.neurips.cc/paper_files/paper/2018/hash/1c336b8080f82bcc2cd2499b4c57261d-Abstract.html (Accessed: 19 November 2023).

2. p(z) is modeled to be a standard Gaussian with diagonal covariance. Would the independence between different z_t allow sudden jumps in the parameters of the system? Would it be better to model it as something like a Gaussian process?

Response: Yes, this independence could allow sudden jumps in the parameters of the system. However, because HyperSINDy learns to approximate the (potentially complex) distribution of states and derivatives (i.e., q_\phi(z|x,\dot{x} is kept close to the prior p_theta(z)), independent samples z_t are transformed into realizations of system coefficients that are similar to those observed in the original data. In other words, over the course of training, HyperSINDy learns to utilize i.i.d. dynamical noise to mimic the dynamical behavior of the observed data – including tendency toward sudden jumps. This is an important feature distinguishing HyperSINDy from other probabilistic methods that do not explicitly model the role of dynamical noise in the generation of the observed time series. We agree that Gaussian processes would provide an ideal framework for modeling the noise, particularly when moving beyond the i.i.d. noise case (see also (Casale et al., 2018; Mishra et al., 2022)). For our work, we assumed that the noise was i.i.d. in time, meaning that the independence of z_ts is a desirable modeling feature (as per the standard SDE framework). However, we note that extension to temporally correlated noise will be an important direction for future inquiry.

3. Related to the last question: does the discretization step size influence the model learning result?

Response: The discretization step size can influence results, although this is not an issue connected to HyperSINDy per se – rather, this relates to the more general and well-studied phenomenon of “finite-time effects” emerging in the numerical analysis of continuous stochastic dynamics based upon discretely sampled time series. Algorithmically correcting for such effects has been the topic of numerous works (e.g., (Ragwitz and Kantz, 2001; Lade, 2009; Honisch, 2011; Boujo and Cadot, 2019)), including one such approach recently introduced for stochastic SINDy (Callaham et al. 2021). Recent works also introduce promising approaches toward finite-time correction in the high-dimensional setting (e.g., (Brückner, Ronceray and Broedersz, 2020; Frishman and Ronceray, 2020)), suggesting a number of promising future directions for HyperSINDy to be extended to the theoretically continuous case.

4. I don't think what "E-SINDy" stands for is ever introduced in the paper.

Response: Thank you for pointing out this oversight! We have updated the Introduction to define E-SINDy.

References

Bauer, S. et al. (2017) ‘Efficient and Flexible Inference for Stochastic Systems’, in Advances in Neural Information Processing Systems. Curran Associates, Inc. Available at: https://papers.nips.cc/paper/2017/hash/e0126439e08ddfbdf4faa952dc910590-Abstract.html

Boujo, E. and Cadot, O. (2019) ‘Stochastic modeling of a freely rotating disk facing a uniform flow’, Journal of Fluids and Structures, 86, pp. 34–43. Available at: https://doi.org/10.1016/j.jfluidstructs.2019.01.019.

Brückner, D.B., Ronceray, P. and Broedersz, C.P. (2020) ‘Inferring the Dynamics of Underdamped Stochastic Systems’, Physical Review Letters, 125(5), p. 058103. Available at: https://doi.org/10.1103/PhysRevLett.125.058103.

Casale, F.P. et al. (2018) ‘Gaussian Process Prior Variational Autoencoders’, in Advances in Neural Information Processing Systems.

Frishman, A. and Ronceray, P. (2020) ‘Learning Force Fields from Stochastic Trajectories’, Physical Review X, 10(2), p. 021009. Available at: https://doi.org/10.1103/PhysRevX.10.021009.

Honisch, C. (2011) ‘Estimation of Kramers-Moyal coefficients at low sampling rates’, Physical Review E, 83(6). Available at: https://doi.org/10.1103/PhysRevE.83.066701.

Lade, S.J. (2009) ‘Finite sampling interval effects in Kramers–Moyal analysis’, Physics Letters A, 373(41), pp. 3705–3709. Available at: https://doi.org/10.1016/j.physleta.2009.08.029.

Mishra, S. et al. (2022) ‘piVAE: a stochastic process prior for Bayesian deep learning with MCMC’, Statistics and Computing, 32(6), p. 96. Available at: https://doi.org/10.1007/s11222-022-10151-w.

Pawlowski, N. et al. (2018) ‘Implicit Weight Uncertainty in Neural Networks’. arXiv. Available at: https://doi.org/10.48550/arXiv.1711.01297.

Ragwitz, M. and Kantz, H. (2001) ‘Indispensable Finite Time Corrections for Fokker-Planck Equations from Time Series Data’, Physical Review Letters, 87(25), p. 254501. Available at: https://doi.org/10.1103/PhysRevLett.87.254501.

Solin, A., Tamir, E. and Verma, P. (2021) ‘Scalable Inference in SDEs by Direct Matching of the Fokker– Planck– Kolmogorov Equation’, in Advances in Neural Information Processing Systems.

Weaknesses

1. Experimental results on higher dimensional datasets might be a bit lacking. One of the important claims of the advantage of HyperSINDy is that it circumvents the curse of dimensionality which hinders the performance of other methods. However, only the HyperSINDy results for one 10D system is given. It might be better if the authors can clarify how the other methods perform on this system, and/or give other examples of high dimensional systems.

Response: Thank you for raising this point. Indeed, the ability to scale to higher dimensions is a distinguishing feature of HyperSINDy. To clarify, the limited scalability of existing SDE discovery methods precludes an assessment of their performance on the 10D system. However, we agree that it would still be informative to benchmark HyperSINDy’s performance against other methods for probabilistic model discovery (of deterministic ODEs) in this high-dimensional setting.

To this end, we have updated the manuscript with new experiments applying E-SINDy to the 10D Lorenz system. In particular, we simulated the 10D Lorenz system with 3 different levels of noise on the forcing term ($\sigma=0$, $\sigma=5$, $\sigma=10$) and trained HyperSINDy and E-SINDy on each trajectory. We evaluated the RMSE of the mean and standard deviation of the discovered coefficients, as well as precision and recall. Refer to Table S5 in the appendix for the full results. Although E-SINDy achieves slightly improved performance in the noise-free case ($\sigma=0$), HyperSINDy outperforms E-SINDy across all metrics on both of the higher noise levels.

Questions

1. In section 3, in "H implements the implicit distribution p(xi_z | z)", why is it the "implicit distribution?" From my understanding, shouldn't xi_z just be a delta distribution (deterministic) on H(z)?

*Response: We apologize for the lack of clarity – we believe this reflects the distinction between z as a random sample vs. a random variable. For a specific sample $\hat{z}, Xi_{z} = H(z)$ is indeed deterministic. For the more general case corresponding to z as a random variable, p(Xi | z) refers to the posterior distribution of SINDy coefficients, conditioned on z. In this setting, the hypernetwork can be interpreted as an implicit distribution for p(Xi | z) (i.e. the neural network weights, biases, and architecture parameterize the distribution p(Xi | z), where z is a random variable). We use “implicit distribution” to mean that we do not actually have full access to the distribution of p(Xi | z), but can still sample from it. For example, if we sample an ensemble of z vectors, then feed that ensemble into H, we are effectively sampling from p(Xi | z). We have revised this line for clarity, including a pointer to (Pawlowski et al., 2018).

Weaknesses

1. Random differential equations are conjugate to stochastic differential equations. It is unclear how to convert a general SDEs into its random ODEs representations, for example, the Langevin type dynamics.

Response: Although this conjugacy is well-established, explicit transformations between SDEs and RDEs are not always straightforward. The standard procedure involves substitution via an Ornstein-Uhlenbeck process, the solution to a linear SDE. Following Han & Kloeden (2017), we may demonstrate this transformation for a scalar SDE with additive noise:

becomes

where $Z_t \coloneqq X_t - O_t$ and $O_t$ is the stationary Ornstein-Uhlenbeck process satisfying the SDE $dO_t = -O_tdt + dW_t$.

A similar conjugation can work for many cases of additive or multiplicative SDEs. In general, though, this conversion is not straightforward to implement. Thus, although we appeal to RDE-SDE conjugacy in order to validate results obtained in the RDE framework, the ability to convert between these representations in practice is more challenging.

On the other hand, this challenge of converting between SDE and RDE formulations highlights unique advantages of our problem formulation. Indeed, much of the interest in RDEs is motivated from the practical advantages afforded by this formulation in certain scenarios in comparison to SDEs (e.g., (Bauer et al., 2017; Solin, Tamir and Verma, 2021)).

We have added this discussion to the Supplementary Appendix.

2. This manuscript lacks a comparison with other methods.

Response: Our original submission included comparisons with the leading “ensemble SINDy” method for robust equation discovery in the presence of noise, and with “stochastic SINDy” for a two-dimensional SDE discovery model. We have updated the manuscript to also include comparisons with Bayesian Spline Learning. Refer to Table S6 in the appendix for comparisons between HyperSINDy, E-SINDy, and Bayesian Spline Learning on a Lotka-Volterra RDE simulation.

3. Although the authors have commented in the manuscript, it is still unclear if this HyperSINDy framework can handle complex noise terms as well as the robustness of noises.

Response: We have performed new experiments to address the question of how our method performs under different scenarios for the noise. Refer to Figure S5 in the appendix for results comparing HyperSINDy and E-SINDy on a Lotka-Volterra system simulated as an RDE, but where each coefficient is drawn from a Half-Normal distribution. In this case, the distribution of coefficients does not match the Gaussian assumption of the latent space; however, HyperSINDy is still able model the complex coefficient distributions due to the expressivity provided by its neural network. This result establishes that the method has potential to generalize beyond the iid noise case, while leaving for subsequent work a more systematic validation in different noise settings.

Questions

1. This manuscript could have been enhanced if it can provide examples of learning underdamped Langevin systems, for example, learn the harmonic oscillator under thermal bath.

Response: Thank you for raising these possibilities for further extensions of our approach. Our updated Appendix clarifies that explicit RDE-SDE transformation is not always trivial, meaning the recasting of these canonical physics SDE problems may not be entirely straightforward; however, we agree that it represents an important future direction for broadening the scope of our framework.

2. In particular, the manuscript could have been enhanced if it can provide an appendix discussion on how to construct a Random ODE representation for a general SDE.

Response: Thank you for making this suggestion, which we agree aids the interpretation of our framework. We have updated the appendix as requested (see also response to Weaknesses #1 above).

3. The manuscript could have been enhanced if it can provide numerical examples when different types of noises are added to the observation data.

Response: Please see response to Weaknesses #3 above.

References

Bauer, S. et al. (2017) ‘Efficient and Flexible Inference for Stochastic Systems’, in Advances in Neural Information Processing Systems. Curran Associates, Inc. Available at: https://papers.nips.cc/paper/2017/hash/e0126439e08ddfbdf4faa952dc910590-Abstract.html

Han, X. and Kloeden, P.E. (2017) Random Ordinary Differential Equations and Their Numerical Solution. Singapore: Springer Singapore (Probability Theory and Stochastic Modelling). Available at: https://doi.org/10.1007/978-981-10-6265-0.

Imkeller, P. and Schmalfuss, B. (2001) ‘The Conjugacy of Stochastic and Random Differential Equations and the Existence of Global Attractors’, Journal of Dynamics and Differential Equations, 13(2), pp. 215–249.

Thank you for your comments on how the manuscript can be made more compelling. We wish to suggest that there are a couple points of confusion here – we apologize for the lack of clarity and have done our best to address these points in the revised manuscript.

Weaknesses

1. It is a typical A+B type of paper. Each part is well studied and author glue them together and demonstrate it in several simple examples. I don't think there is enough novelty here.

Response: We believe the present contribution goes well beyond a “typical A+B paper”, in the sense that the two methods claimed to be simply glued together (hypernetworks (Ha, Dai and Le, 2016; Pawlowski et al., 2018) and SINDy (Brunton, Proctor and Kutz, 2016; Kaptanoglu et al., 2022) are just two pieces of an overarching conceptual and technical framework being put forward. Moreover, these two methods refer to high-level frameworks/architectures with many possible implementations (e.g., a hypernetwork is simply a neural network that predicts the weights of another network – a concept we reimagine for predicting the coefficients of an analytical expression). Thus, the integration of the two methods (which have arisen in different academic communities for very different purposes) is far from straightforward.

At a conceptual level, integration of these methods is entirely contingent on our innovative problem formulation in terms of random differential equations. This alternative formulation for stochastic dynamics enabled us to leverage the unique advantages of numerous algorithms/techniques (including SINDy and hypernetworks, but also, e.g., variational autoencoders and backpropagation-compatible L0 loss) toward our overarching goal of stochastic dynamical modeling.

At a technical level, unifying all these approaches within a common framework is a highly non-trivial task, requiring strategies to appropriately balance multiple losses relating to equation sparsity, overall model complexity, reconstruction error, and proximity to a Gaussian prior (that is, promoting discovery of a model approximating a canonical SDE).

All in all, while it can be heuristically useful to summarize the proposed approach as combining hypernetworks and SINDy, this intuitive picture belies a much more sophisticated integration of techniques and ideas. We believe that this claim is consistent with evaluations from the other reviews, who each note the novelty/innovation of the proposed framework.

2. All three examples are artificially made for this algorithm. All examples are corrected identified but I am not impressed unless authors are able to demonstrate some non-trivial RODE. The second example equation (11) is not even a valid example of stochastic Lotka-Volterra. I don't know what is N(0,1) on the Right hand side means here.

Response: We believe that this impression stems in part from notation differences between fields. We have modeled systems subject to white noise (thus, a general setting for stochastic dynamics), though we adopt notation more familiar in statistics/machine learning. Thus, N(0,1) denotes a Gaussian distribution with mean=0, standard deviation=1; iteratively sampling from this distribution amounts to simulating a specific realization of a driving white noise process (which is defined by random variables distributed according to N(0,1) at each time point t (Han and Kloeden, 2017). As such, our selected problems are consistent with widely used stochastic modeling problems (namely, classic SDEs) in which the source of stochasticity is taken to be a Gaussian white noise. Nonetheless, we agree that this notation can obscure connections to standard SDE representations; as such, we have rewritten equation 11 in canonical SDE matrix notation involving a 2D Wiener process (which we simulate via the Euler-Maruyama integration scheme), thus clarifying that it is indeed a valid Lotka-Volterra SDE.

3. Authors have limited knowledge on RODE here in fact not all SDE can be transformed to RODE and vice versa. And in general RODE case, z is not independent with x.

Response: As we understand, any finite-dimensional, ordinary SDE can be transformed to an RDE; we agree that the reciprocal is more restricted to RDEs driven by stochastic processes that are the solutions of SDEs (Han and Kloeden, 2017). To our knowledge, this reciprocity is well-established since the results of (Imkeller and Schmalfuss, 2001), and has been usefully invoked in other numerical investigations into stochastic dynamics where an RDE formulation carries practical advantages (e.g., (Bauer et al., 2017)). We agree that in the general case, the noise process z(t) need not be state-independent. However, the state-independent case (i.e., an autonomous noise process) represents the most common problem setting for stochastic dynamics, and thus carries broad relevance beyond the specific models studied herein.