For better readability, we restate your comments/questions in italic, our response(s) are marked by ▶

This work resonates closely with inverse problems for dynamical systems, where several works have been conducted to infer the initial parameters ...

▶ Currently, we discuss related work on the inverse problem, but only in the context of models of collective behavior. We thank the reviewer for suggesting a broader perspective and comparing with prior learning-based works on inverse problems for dynamical systems in general, including the three listed ones. This will certainly be a good addition to the manuscript.

Regarding refs 1-3, we note that the underlying differential equations (i.e., classic PDEs) differ from the multi-particle systems we consider. In our case, forces on each particle depend on its relative position towards all other (or its neighboring) particles, resulting in strong couplings between individuals. For specific systems (Volume exclusion), it is possible to derive parametric PDEs that describe the asymptotics of infinitely many particles (Oelschläger, 1990), and ideas from references 1-3 might be applicable. However, in general, the number $N$ of particles influences the dynamics. E.g., for the D'Orsogna model, particle distances can converge to 0 as $N\to \infty$.

I think the paper should include some additional discussion about the computational limitations of the method ...

▶ We agree; please see our General Response section.

Despite the capability of their model to incorporate persistent homology, enabling the use of topological features, they only utilize classical persistent homology ...

▶ Vietoris-Rips persistent homology is the de-facto workhorse for point cloud data. While one could use other geometric complexes (e.g., Alpha complexes as mentioned by reviewer dbj5) or different vectorization strategies (ATOL, persistence images, etc.), we decided to work with the most prevalent tools, but acknowledge that other choices are possible. Also, using Vietoris-Rips PH facilitates straightforward comparisons to previous topology-related works on the same problem, such as Crocker Stacks or the PSK approach.

Importantly, our goal was not to assess the impact of different design choices but to introduce a latent dynamical model based on topological per-time-point descriptors for the inverse problem at hand. We are convinced that future work will explore different variants of our approach that may be more suitable in certain situations (e.g., using LSTMs in case of discrete dynamical systems, see 83FG).

Does the model output a single point estimate for the parameters of interest or a corresponding distribution (mean and std)?

▶ At the moment, all models output point estimates of the sought-for parameters. However, one could easily sample $q_{\\theta}$, integrate forward, and obtain a distribution of estimates. We will add tables to the appendix that report an estimated mean & std. dev. from this sampling.

Can authors describe the regression objective w.r.t to underlying parameters of interest, ...?

▶ The objective is to predict the simulation parameters that led to a particular realization of a point cloud sequence. Take the PH-only variant v1 of Fig. 1, for instance: input is a sequence of vectorized persistence diagrams $v_{\\tau_0}, \\ldots, v_{\\tau_n}$. The encoder Enc_\\theta yields parameters of the approx. posterior q_\\theta from which we sample and integrate the latent ODE forward in time to get latent states along the trajectory; this latent state sequence is then fed through Enc_\\alpha which summarizes the sequences into a vector and linearly maps the latter to simulation parameter estimates $\\hat\\beta = (\hat{\beta}\_1, \ldots, \hat\beta\_P)$. During training, we minimize the mean-squared-error (MSE) between the predictions and the ground-truth simulation parameters, implicitly making a Gaussian noise assumption. One could write the regression model as $\\hat\\beta = \\phi(v_{\\tau_0}, \\ldots, v_{\\tau_0}) + \epsilon, \epsilon \sim \\mathcal N(0,\sigma I)$ where $\phi$ subsumes all steps above.

Since PH does not account for the temporal changes in the point cloud, how is the work different from this https://arxiv.org/abs/2406.03164? ...

▶ As correctly pointed out, we do not track topological features on the level of individual points in a persistence diagram but on the level of vectorized persistence diagrams. The crucial difference to the arXiv paper is that the integration of neural ODEs is on the message passing level of a GNN for a fixed graph. In our context, this means that the graph object needs to be such that each vertex corresponds to the same particle at all times. However, such correspondences are unknown or hard to obtain in practice. In fact, in our work, we deliberately avoid this "tracking" step.

Does the fitting of the ODE system deal with stiffness as you are fitting the topological features via an ODE, which can be non-smooth in between, leading to inefficient modeling?

▶ This is an interesting point! However, we do not model the dynamics directly in the space of persistence diagrams but rather in a latent space learned from diagram vectorizations. This modeling could be inefficient to some extent, yes, but the latent ODE approach seeks to learn the most suitable latent space for the task at hand, which is to minimize prediction error for the model parameters (via the MSE) and to minimize the reconstruction error for the vectorized persistence diagrams (as part of the ELBO). This strategy is common in the literature, as, e.g., seen in the PhysioNet 2012 experiments of (Rubanova et al., 2019), where variables include binary indicators (e.g., of whether mechanical ventilation is used at available time points).

For better readability, we restate your comments/questions in italic, our response(s) are marked by ▶

I think that the sentence in line 32 "... due to missing correspondences between individuals..." is hard to read. Which missing correspondences?

▶ By "correspondences" we mean that many previous works rely on being able to unambiguously identify particles across point clouds over time. We will clarify this unfortunate formulation.

The proposed approach does not track of the evolution of individual topological features of the persistence diagram. This could be important in scenarios requiring fine-grained topological information.

▶ Indeed, our method does not retain information about how individual homology classes evolve. We model the evolution of features on the level of persistence diagram vectorizations, not on the level of individual points in the persistence diagram(s). Persistence vineyards would allow us to do this, but they are (i) hard to vectorize for ML tasks, (ii) computationally expensive to compute, and (iii) often infeasible to obtain without particle correspondences over time. Following your suggestion, we will extend our discussion to include these points.

The authors claim that their method is more efficient than other methods, but computation times are not reported in the main text.

▶ Please see our runtime analysis in the General Response section.

In the persistent homology paragraph (starting at line 154), the authors state that there is a decomposition in persistence diagrams (births and deaths) when using abelian groups ...

▶ When sketching persistent homology in Sec. 3, our writing unintentionally and wrongly suggests that birth/death times (interval compositions) are well-defined for homology with coefficients in arbitrary Abelian groups, e.g., in the presence of torsion. We will clarify this part; thank you for pointing it out. We use coefficients in the field $\\mathbb{Z}/2\\mathbb{Z}$ (as is usual).

Only one vectorization method is used, without significant justification. I would have appreciated a more detailed ablation on this point.

▶ Our contribution is not within the particular realization of our core idea, but in the surprising insight that one can obtain remarkable predictive performance from capturing the dynamics of vectorized topological summaries per time point, without tracking particles, or without computationally-expensive summaries of time-varying persistence. On could, however, use other vectorization techniques; we experimented with ATOL vectorizations (which are conceptually very close to ours) and persistence images, but only observed minor changes in predictive performance. Nevertheless, ATOL vectorizations are not stable, and persistence images can be hard to parametrize (and discretize) without obtaining (very) high-dimensional descriptors.

Maybe I'm wrong, but when I access the arxiv version of reference [23] in your paper, Table 4 contains execution times, and not SOTA results.

▶ Thank you for pointing this out! Table 4 is meant to refer to the table in the submitted manuscript, not in the referenced work. We will rephrase our formulation.

Why do you use 𝑅2 ? I thought that it was a misleading metric [4]. This is my main concern, and the core reason why the score is a borderline accept.

▶ We agree, $R^2$ has known flaws, much like other metrics to assess regression performance, yet it is easy to interpret. As a summary impression across all predicted parameters, MSE and RMSE cannot be applied due to differing ranges of simulation parameters. We used $R^2$ and SMAPE to have two different evaluation scores. Our joint variant (v3) outperforms the state-of-the-art in both of these scores across all datasets, and even our PH-only variant (v1) does so on 3 out of 4 datasets.

To further reduce evaluation bias, we will add the explained variance (EV) as suggested in your ref [4] to all tables, and also list the RMSE per parameter in our appendix. Below is an example for the dorsogna-1k data (parameters $C,l$ are as in (Guisti et al., 2023)):

Can your approach be used for discrete dynamical systems (where time is discrete and not continuous?)

▶ We are confident that our method can be used with discrete systems. However, some components of our architecture may need to be adjusted; e.g., switching from neural ODEs for modeling latent dynamics to, say, LSTMs or GRUs. Thank you for this comment, we will add a remark!

Since the method is very general and can work with any kind of vectorizations, not just topological ones, is there a way to regularize or add explicit bias during training ....?

▶ If we understand your question correctly, you are asking whether some form of topological prior could be used to inform the neural ODE training process. This is an interesting question, and while we have not experimented with topological regularization so far, it might be possible to enforce specific topological/geometrical properties through an appropriate loss (on persistence diagrams) and differentiating through the PH computation (as done in prior works).

For better readability, we restate your comments/questions in italic, our response(s) are marked by ▶

ad Clarity: We tried to balance the presentation of our conceptual idea vs. the technical detail of its particular realization. We can provide more detail on the ELBO objective and the added regression loss in the appendix. The same holds for the PD vectorization technique we use. As both parts appeared separately in prior work, we sacrificed technical detail for those parts in favor of a more comprehensive conceptual perspective. We will add these details to the appendix to facilitate reproducibility, as you suggested, including pointers to the respective parts in the source code.

Further, for training, the regression loss is the MSE between the parameter estimates and the ground truth simulation parameters. One may also use $R^2$ (or any other reasonable objective) during training, but we found the MSE to already yield good results. During inference/evaluation (on held-out testing data), we report the $R^2$ and SMAPE computed from predictions vs. ground truth values, just as you correctly inferred.

Related to my comments above, could you make explicit the way the parameters of the systems are encoded and learned? Say that I consider the Volume exclusion model with parameters. How are these parameters actually estimated once Enc_\\alpha, Enc_\\theta, Dec_\\gamma are trained with the ELBO + regression rule. How are the $\hat{\beta}_t$ (showcased in Figure 3) used afterward?

▶ Inferring the simulation model parameters (e.g., for the Volume Exclusion model) is done as follows: For a given sequence of point clouds (at time points $\tau_0, \\dots,\\tau_n$), one first pre-computes Vietoris-Rips persistent homology (e.g., for $H_0, H_1, H_2$) per available time point and vectorizes the diagrams. This yields a sequence of vectors, i.e., one vector per available time point. This sequence is then fed to the encoder $Enc_\theta$ which outputs a parametrization for the approximate posterior q_\\theta(z_{t_0}|\\{ v_{\tau_i} \\}). Upon sampling from q_\\theta, we get an initial latent state $z_{t_0}$ and integrate the latent ODE forward in time to $t_n$. Enc_\\alpha then summarizes this sequence of latent states and linearly maps it to a vector of simulation parameter estimates $\beta_1, \\ldots, \beta_P$. (For Volume Exclusion, we have four simulation parameters, so $P=4$).

(more genuine question out of curiosity) In some sense, the seminal point cloud contains implicitly all the possible topological information one could extract using persistent homology. It also seems that PointNet-like architectures ...

▶ An important point in this context is that RipsNet is trained (in a supervised manner) to predict precomputed vectorizations of persistence diagrams, not the persistence diagrams themselves. This means the desired vectorized diagram is known in advance and can be used in a loss function. Clearly, this strategy is designed to guide the network towards capturing topological information. Instead, in our experiments, the PointNet++ directly predicts the parameterization of the simulation; hence, it would have to learn topological information implicitly (as there is no external guidance towards this goal). The only guidance given is a loss wrt. the predicted simulation parameters. It might, however, be possible to replace our full pre-computation step (for Vietoris-Rips PH) with a pre-trained RipsNet (i.e., pre-trained on point clouds from such simulation experiments). We did not explore this direction here, but we greatly appreciate the comment.

Still related to the above, the paper mentions that it uses PD with homology dimension 0, 1 and 2. Did you run the experiments using $H_0$ alone? If the performance significantly decreases, that would be a nice example where non-trivial topology is useful (as frustrating as it can be, I often observed that $H_0$ alone, i.e. looking mostly at the distribution ...

▶ Yes, initially, we experimented with $H_0$-only and found that prediction performance drops in that case. The situation is less clear when also including $H_2$ features, where the results are more mixed and often not noticeably different. From a more quantitative perspective, below you can find a table for the dorsogna-1k experiment comparing $H_0$ vs. $(H_0, H_1)$ vs. $(H_0, H_1, H_2)$ for our approach (i.e., the v1 variant from Fig. 2):

In the conclusion, you mention that computing $H_2$-PD with Rips filtration is quickly computationally expensive (I completely agree), and propose to use alpha-filtration instead. Is there any difficulty to do so right now? (Alpha filtrations are readily implemented in several Tda-libraries according to https://cat-list.github.io/ ).

▶ There is no inherent limitation to switching out Vietoris-Rips complexes for Alpha complexes. However, we did find in a preliminary experiment that different implementations of Alpha complexes (some of which are in the link you provided) surprisingly yield different diagrams, and without any further in-depth investigation refrained from using the latter in our experiments. Nevertheless, switching from VR to Alpha complexes would, as you suggested, significantly improve the performance of the preprocessing step in which we compute PH for different dimensions.

For better readability, we restate your comments/questions in italic, our response(s) are marked by ▶

Ad contribution: The reviewer is correct in that we use the latent ODE framework of (Rubanova et al., 2019). However, the latter is only one particular variant (of many; other options are, e.g., latent SDEs or Continuous Recurrent Units (CRUs)) for modeling latent dynamics. Our main contribution is the idea of modeling the point cloud dynamics through the lens of vectorized topological summaries, i.e., capturing topological changes on a persistence diagram level, as opposed to capturing changes on the level of individual points in such topological summaries. Even this arguably "simple" approach yields remarkable predictive performance for the parameters of governing equations of collective behavior. Furthermore, our approach avoids the need to track particles over time (and infer velocities) or to rely on computationally expensive summaries of time-dependent persistence.

Ad empirical data: Thank you for that comment. Yes, our approach would work on the data from (Bialek et al., 2012). In fact, one motivation for our work was to alleviate the challenge of having to compute/estimate per-particle velocities (this is required in the work of Bialek et al. as they compute correlations from normalized velocities). In many real-world settings, one might not even be able to unambiguously track particles (let alone when the cardinality of the point clouds may vary as well). Instead, our approach only hinges on the position of the particles (and can even naturally deal with point clouds of varying size, as demonstrated with the volex-10k experiment).

Ad training time: Please see our runtime analysis in the General Response section.

Neural ODEs are well known to be very difficult to train (see e.g. Dupont et al. (2019), Finlay et al. (2020), Choromanski et al. (2020) or Pal et al. (2021), just to cite a few). Can you comment on how difficult (or easy) was to train your Latent ODE network on the synthetic data you studied? What about the training time?

▶ First, we did not experience any difficulty training the full model. However, we did only experiment with the arguably simplest ODE model (an autonomous ODE). We did preliminary experiments also with a non-autonomous variant, but did not observe any noticeable improvements. Nevertheless, our experiments support the conclusion that even the "simplest" choice of latent ODE already suffices to largely outperform existing techniques.

How did you parametrize the decoder network? I don’t find it in the manuscript.

▶ We apologize for not being clear enough on that point. The decoder (aka reconstruction network, denoted as Dec_\\gamma in Fig. 2 of the manuscript) is a simple 2-layer MLP with ReLU activation that maps latent states (from the latent ODE) to reconstructed persistence diagram vectorizations. We will include this information in the revised manuscript's appendix (and point to the locations in our source code).

Did you consider adding some noise to your observations in point cloud space? As presented, your network should easily be able to handle noisy observations. Similarly, how does the persistent homology representation deal with noise? These questions are of course relevant if one wants to apply your methods to empirical data.

▶ Currently, only the governing equations of the Vicsek model (see Fig. 3) include noise through the Brownian motion. In general, we point out that due to the stability of the persistence diagrams wrt. perturbations of the input (see the Stability/Continuity aspects section of the manuscript), our method is suitable in case of observation noise (unless the noise is excessively large). E.g., on dorsogna-1k where points are within $[-1,1]^3$, adding Gaussian noise $\\mathcal{N}(0,\\sigma)$ to all point coordinates and all time points, ($R^2$,SMAPE) drops from (0.846 $\\pm$ 0.011, 0.097 $\\pm$ 0.005) to (0.822 $\\pm$ 0.016, 0.098 $\\pm$ 0.002) at $\\sigma=0.01$, and to (0.730 $\\pm$ 0.013, 0.150 $\\pm$ 0.001) at $\\sigma=0.1$.

How many points from the latent path are used as input to the regressor model? Does the model work with a single point, as e.g. the last point along the latent trajectory? It’d have been nice to understand what information of the latent path is important for the regression task, specially given that completely dispensing from the latent dynamics still gives compelling results (i.e. Table 3).

▶ Again, we apologize for not being precise enough. In fact, we experimented with multiple variants: you are correct that one could use the “last” point along the latent trajectory. Another variant would be summarizing the latent trajectory via a signature approach (although this summary can become quite high-dimensional). The most conceptually elegant approach, from our perspective, was to re-use the encoder architecture (i.e., a duplicate of the mTAN encoder module with its own set of parameters, denoted as $Enc_{\alpha}$ in Fig. 2). However, instead of receiving an unequally-spaced sequence of vectorized persistence diagrams, $Enc_{\alpha}$ receives ALL points along the latent trajectory. In our case, “all” means 100 latent states at equally spaced time points in $[0,1]$, which we obtain by integrating the latent ODE forward in time.

Also, our baseline (w/o dynamics) could be considered quite strong since it uses an mTAN encoder (i.e., an attention-based approach) that directly maps sequences of vectorized persistence diagrams to parameter predictions. While mTANs can attend to relevant information in a sequence, this baseline still consistently falls short of our approach, which explicitly models the dynamics.

I thank the authors for their detailed responses. After reading the reviews from other reviewers and your replies to them, I have decided to increase my score.

I believe that, if the updated manuscript addresses all the points raised during this discussion session, it will make a valuable contribution.