No Equations Needed: Learning System Dynamics Without Relying on Closed-Form ODEs

(C) Bifurcations

We believe that our framework is uniquely positioned to perform quite well on systems exhibiting bifurcations (when a small change to the parameter value causes a sudden qualitative change in the system's behavior). In our framework, bifurcation occurs when the composition map predicts a different composition. As discussed (now more thoroughly) in Appendix G.1, in the future, the semantic predictor may take as input not only the initial conditions but also other auxiliary parameters. We can then represent the composition map as a decision tree that divides the input space into different compositions. This decision tree then informs us where bifurcations occur.

We hope the following proof of concept based on the current implementation demonstrates that it is a viable approach. Instead of predicting a trajectory from its initial condition, we fix the initial condition to be always the same and predict a trajectory based on the parameter $r$ that we observe in our dataset. We generate the trajectories given the following differential equation $\dot{x} = rx - x^2$ the initial condition $x(0)=1$ and $r$ sampled uniformly from $(-1,2)$. We choose the set of compositions to be $(s_{+-h}),(s_{-+h})$ and record the position of the bifurcation point found by our algorithm (as opposed to the ground truth). The mean absolute error for different noise settings can be seen in Appendix B.9 and the table below. Note that the range of values of the trajectory is $(0,2)$, and even in high noise settings, the location of the bifurcation point can be identified.

---

Questions

Q1: Please see our response in "Interpretability and generalizability"

Q2: Please see our response in "Experiments on real datasets".

Q3: Please see our response in "Bifurcations".

Q4: Partial differential equations (PDEs) no longer describe trajectories (a function defined on the real line) but rather fields that are defined on multi-dimensional surfaces such as planes. Extending direct semantic modeling to PDEs requires developing a new definition of semantic representation of a field. Moreover, as the initial condition is no longer a single number, the semantic predictor needs to take whole functions as inputs (for both initial and boundary conditions). We have added this discussion to Appendix G.4.

Q5: We have now updated all differential operator symbols to upright type.

Dear Reviewer BR8f,

Thank you very much for your review. We are glad that you found our approach to be novel. By addressing your questions and concerns, we truly believe we have strengthened our paper. We provide our responses below and point to relevant changes in the paper.

---

(A) Interpretability and generalizability

Interpretability

The main advantage of our approach in terms of interpretability is the fact that we do not have to perform any mathematical analysis to obtain a semantic representation of the model (as it is defined in the paper). Obtaining the same information from a closed-form ODE is often possible (if the equation is not too complex) but may be time-consuming or require mathematical expertise.

Generalizability

Closed-form ODEs may generalize very well if the discovered equation is the correct one. However, finding the correct equation may be hard or sometimes even impossible if the phenomenon is not governed by an ODE. For instance, the dynamics may depend on history (e.g., delay and integral differential equations), involve unobserved variables (e.g., pharmacological model), or empirically determined functions (e.g., stress-strain curve). In Section 6.2, we showed how ODEs may fail to generalize outside of the training time domain and how, by editing and incorporating semantic inductive biases, we can increase the generalizability of Semantic ODE to the unseen time domain.

We have included another experiment in Appendix B.10, showcasing how we can generalize the model in Section 6.2 to initial conditions $(x_0)$ not seen during training. We observe that each property function in Figure 8 looks approximately like a linear function. Thus, we fit a linear function to each of these functions and then evaluate our model on initial conditions from range $(1.0, 1.5)$. Note that our training set only contained initial conditions from $(0, 1)$. We compare the performance of this model to ODEs from Table 2. The results can be seen below. Our model has suffered only a small drop in performance even though it has never seen a single sample from that distribution. It also performs much better than any other ODE tested.

---

(B) Experiments on real datasets

Thank you for your suggestion regarding additional experiments. We have now added experiments on two real datasets (on tumor growth and drug concentration). We added these results to Appendix B.8 and reproduced them below. "Semantic ODE*" is a variant of Semantic ODE where we incorporated a semantic inductive bias about the shape of the trajectory. We specified the composition to always be $(s_{-+c},s_{++u})$ for the tumor growth dataset and $(s_{+-c},s_{--c},s_{-+h})$ for the drug concentration dataset. Note that the implementation of WSINDy we used cannot work with trajectories as sparse as the drug concentration dataset.

(D) Minor issues

Thank you for spotting the typos. We have fixed them now.

---

Questions

Q1: Thank you for asking this question! We should have mentioned this in the implementation details. Indeed, in Algorithm 1, before we decide to switch from one composition to the other, we verify that both the number of samples and the interval length are above user-defined thresholds. In experiments, we require at least two samples and the interval length to be at least 10% of the size of the entire domain of the map. We have now included this in Appendix C.1. There is also a maximum number of splits decided by the user. In all experiments, we allow for at most 3 subsets/property maps. The property maps themselves are described as a linear combination of basis functions (each described by just a few parameters), so, in principle, they require fewer samples than models with a larger number of parameters.

Q2: No, each $F_{\text{prop}}$ is trained separately for each composition. They do not share any information. However, we believe some amount of sharing may improve performance where the two adjacent compositions are similar, and the exact boundary is uncertain. We consider this an interesting extension for future work. We have included it in Appendix G.6.

Q3: Before predicting with $F_{\text{traj}}^2$, we ensure that the predicted trajectory's transition points are at most $\epsilon$ from those dictated by the property map. If we cannot find such a trajectory, then we default to using $F_{\text{traj}}^0$. In our experiments, we choose $\epsilon=0.001$. In general, we consider $F_{\text{traj}}^0$ to be an approximation of $F_{\text{traj}}^2$ used for backpropagation. We also have a penalty term added to the training loss that encourages smoother trajectories so that finding a $\mathcal{C}^2$ trajectory is easier. We use $F_{\text{traj}}^2$ for both validation and test errors, so we do not expect any major issues associated with that.

References

[1] Chen, R. T., Rubanova, Y., Bettencourt, J., & Duvenaud, D. K. (2018). Neural ordinary differential equations.

[2] Holt, S. I., Qian, Z., & van der Schaar, M. (2022). Neural laplace: Learning diverse classes of differential equations in the laplace domain.

[3] Lu, L., Jin, P., & Karniadakis, G. E. (2019). Deeponet: Learning nonlinear operators for identifying differential equations based on the universal approximation theorem of operators.

Dear Reviewer DbZG,

Thank you very much for your review. We are glad you found our work to be completely new and different from previous approaches. By addressing your questions and concerns, we truly believe we have strengthened our paper. We provide our responses below and point to relevant changes in the paper.

---

(A) Capitalization of the title

It seems to us that the current capitalization follows common capitalization guidelines, including AP Stylebook and APA Style. Please let us know if you believe we should change the capitalization of some words to make it consistent.

---

(B) Comparison with black-box approaches

We have now added a comparison with three black box approaches: NeuralODE [1], NeuralLaplace [2], and DeepONet [3]. All results are now included in Table 3 and can also be seen below. In nearly all settings, Semantic ODE outperforms the black-box method. This is likely because we operate in a low-data regime, where we always have fewer than 200 samples measured at most 20 times throughout the trajectory.

---

(C) Downstream tasks

The main difference between direct semantic modeling and traditional two-step modeling is the fact that the user can interact directly with the semantic representation of the model. That not only avoids (often demanding) the analysis process but, more importantly, allows for incorporating semantic inductive biases (Section 6.1) and model editing (Section 6.3), which can significantly increase its extrapolation properties to the unseen time domain (as shown in Section 6.3).

We have included another experiment in Appendix B.10, showcasing how we can generalize the model in Section 6.2 to initial conditions $(x_0)$ not seen during training. We observe that each property function in Figure 8 looks approximately like a linear function. Thus, we fit a linear function to each of these functions and then evaluate our model on initial conditions from range $(1.0, 1.5)$. Note that our training set only contained initial conditions from $(0, 1)$. We compare the performance of this model to ODEs from Table 2. The results can be seen below. We can see that our model has suffered only a small drop in performance even though it has never seen a single sample from that distribution. It also performs much better than any other ODE tested.

Dear Reviewer DbZG,

Thank you once again for your thoughtful feedback, which has significantly improved our paper. We have carefully addressed your comments in our rebuttal and hope you find our responses satisfactory.

We wanted to check if you have any remaining concerns or questions we can address. If our revisions meet your expectations, we hope you might consider increasing your score.

We truly appreciate your time and invaluable insights!

Kind regards,

Authors

Computational cost

We have now included the following discussion in Appendix G.2.

The computational cost of extending semantic modeling to $M$ dimensions depends on the exact approach taken. We will try to make our best guess, assuming the approach remains as similar as possible to the one we describe in the paper.

With our current ideas, the added cost of having multiple dimensions comes mostly from the need to perform the same tasks $M$ times. For $M$ trajectories, we need to fit $M$ composition maps and $M$ property maps.

We do not expect the time to fit a single composition map to increase significantly. The main computational burden of the current implementation comes from fitting a composition to each sample to see how well it fits. For the same number of samples, the computation cost is going to be the same. Although there will be an added cost to find optimal decision boundaries in multiple dimensions (as opposed to one), we do not think this will be a significant burden given current efficient classification algorithms.

Fitting of each property map will require $M$ times the number of parameters. However, it is important to note that the current property maps have just tens of parameters in most cases, so this number is still likely to be a few hundred at most. The number of property maps needed to be fitted (for each of the compositions) may be larger. However, a user would still like to narrow it down to a reasonable number (likely less than 10) to make sure that the composition map is understandable.

Overall, we suspect the training time to increase at least linearly with the number of dimensions, but we do not expect any combinatorial explosion as often happens in symbolic regression, so for the number of trajectories, we may be interested (probably less than 10), the computational costs should not be a major concern.

---

(B) Critical points

As we mentioned in the limitations, our approach is limited to trajectories with a finite composition. We have now elaborated in Appendix G.1 how this can be extended to oscillatory/periodic trajectories by meta-motifs and their associated meta-properties. The idea is to take the infinite composition (infinite sequence of motifs) and represent it using some finite representation. For instance, a trajectory $\sin(t)$ would be described as $(s_{+-c},s_{--c},s_{-+c},s_{++c})$ repeating forever, where the meta-properties may include the "frequency" and "amplitude" that may also vary with time and be itself described using compositions. With this extension, we do not expect the performance to be limited for short-wavelength features as long as the measurements of the trajectory are not too coarse.

However, we admit that the current implementation may struggle if the ground-truth composition is finite but long. In that case, it may fall outside of the chosen set of compositions. We have now explicitly stated this limitation in Appendix G.1.

---

(C) Writing

Thank you for your comments on how to improve our presentation. As suggested, we have now moved Fig. 10 to the main text and expanded the caption of Figure 1.

Dear Reviewer DLBF,

Thank you so much for your time spent reviewing our paper. By addressing your questions and concerns, we truly believe we have strengthened our paper. We provide our responses below and point to relevant changes in the paper.

---

(A) Higher dimensions

Extension

We have now elaborated on extending semantic modeling to multiple dimensions in Appendix G.2. The easiest way to model multiple trajectories (or multi-dimensional trajectories) is to model each trajectory dimension independently but conditioned on the multi-dimensional initial condition $x_0\in\mathbb{R}^M$, i.e., for $M$ trajectories, we want to fit $M$ forecasting models, $F_1, \ldots, F_M$ such that for each $m\in[M]$, $F_{m}:\mathbb{R}^M \to \mathcal{C}^2$. Each $F_m$ can still be represented as $F_{\text{traj}} \circ F_{\text{sem}}^{(m)}$, where $F_{\text{traj}}$ is the same as described in the paper and $F_{sem}^{m} : \mathbb{R}^M \to \mathcal{C}\times\mathcal{P}$. $F_{sem}^{m}$ can similarly be represented using a composition map and property maps for each predicted composition. However, these maps are no longer described by univariate functions. We currently envision the composition maps to resemble decision trees (or similar structures) and the property maps as generalized additive models to make the whole structure comprehensible. Moreover, the semantic predictors can also accommodate auxiliary variables beyond the initial conditions. For instance, the parameters of the ODEs. This approach should scale to multiple dimensions as it is very modular. Each $F_{\text{sem}}^m$ can be analyzed independently. Each property map of $F_{\text{sem}}^m$ corresponding to a particular composition can be analyzed independently. And finally, each shape function for each of the predicted properties can be analyzed independently. We believe this modularity is one of the reasons why even systems with multiple dimensions can remain understandable and, more importantly, can be edited as each change has a localized impact. In contrast, changing a single parameter in a system of ODEs may change all the trajectories in ways that may be difficult to predict without a careful analysis.

Proof of concept

This extension deserves a separate paper, but we wanted to demonstrate a proof of concept. We implemented the property maps as we defined above (each property described as a generalized additive model) and fitted data following a SIR epidemiological model for different initial conditions. To simplify the problem, we prespecify the composition map (we only train property maps). We assume $S$ follows $(s_{--c},s_{-+h})$, $I$ follows $(s_{++c},s_{+-c},s_{--c},s_{-+h})$, and $R$ follows $(s_{++c},s_{+-h})$. We have included the results in Appendix B.5, which shows the predicted trajectories and some of the property maps. The average RMSE on the test dataset is 0.019 for $S$ trajectory, 0.011 for $I$ and 0.014 for $R$. Note that the irreducible error on this dataset (caused by added Gaussian noise) is 0.01. The shown shape functions let us draw the following insights about the model:

• The time when $I$ is at its maximum $t_{\text{max}}(I)$ is on average just below $0.2$. $I_0$ has a relatively large impact on $t_{\text{max}}(I)$ by increasing it by 0.1 for very low $I_0$ or decreasing it by 0.05 for very high $I_0$. The larger the $I_0$, the faster the maximum is achieved.
• $S_0$ also has a negative impact on $t_{\text{max}}(I)$ but it is much smaller ($\pm 0.02$).
• The maximum of $I$ (denoted $I_{\text{max}}$) increases linearly with both $S_0$ and $I_0$. This time $S_0$ has slightly bigger impact ($\pm 0.1$) compared to $I_0$ ($\pm 0.04$) .
• In both $t_{\text{max}}(I)$ and $I_{\text{max}}$, the impact of $R_0$ is insignificant.
• The horizontal asymptote of $R$ increases linearly with all three initial conditions. In particular, the shape function associated with $R_0$ has a unit slope as expected.

An interesting advantage of our approach is that even though three variables describe the system, we do not need to observe all of them to fit the trajectory (similarly to the pharmacokinetic example in the paper). ODE discovery methods assume that all variables are observed, which constrains their applicability in many settings.

Dear Reviewer DLBF,

Thank you once again for your thoughtful feedback, which has significantly improved our paper. We have carefully addressed your comments in our rebuttal and hope you find our responses satisfactory.

We wanted to check if you have any remaining concerns or questions we can address. If our revisions meet your expectations, we hope you might consider increasing your score.

We truly appreciate your time and invaluable insights!

Kind regards,

Authors

(E) Granularity of information

Thank you for drawing our attention to this very interesting point about the granularity of information. We have now included the following discussion in Appendix G.3

We believe the granularity of information (in terms of the shape of the trajectory, its maxima, minima, and asymptotic behavior) should be comparable between direct semantic modeling and other approaches, as long as all approaches are framed as solutions to finding a forecasting model. The key difference is that direct semantic modeling gives immediate access to the semantic representation of the model without any further analysis and allows us to easily edit the model and incorporate semantic inductive biases during training.

There is, however, one way in which ODEs may offer more granular information by assuming a particular causal structure. Note: it may not be correct, and it makes them less flexible. For a discretized ODE, we can write the underlying causal graph as $x(t_0) \to x(t_0 + \Delta t) \to \ldots \to x(t)$. At the same time, our model assumes a more general $x(t_0) \to x(t)$. This generality is useful (as we show in Section 6.4), and it does not negatively impact the performance. However, this model does not allow us to do interventions where we change the value of $x$ at a particular $t$ and predict the future. ODEs do allow for this kind of inference, but of course, there is no guarantee that it is correct. For instance, a system governed by a delayed differential equation has a different causal structure. We believe direct semantic modeling in the future can be extended to accommodate some kinds of intervention, where they can be described as additional features beyond the initial conditions.

---

(F) Chaotic behavior

We could not find chaotic systems that would not have some kind of oscillatory behavior (i.e., that could be described by a finite composition). As such, our current implementation of Semantic ODEs would not be able to successfully predict the behavior beyond the seen time domain. However, we could use it to predict on a bounded time domain. That is why, following your suggestion, we have included an experiment with the Duffing oscillator. Results are shown below and in Appendix B.7. Although our method achieves decent results, this performance would drop drastically if we tried to predict beyond the training time domain. We have clarified this limitation in Appendix G.1.

In general, Semantic ODEs can model significant changes in the trajectory following a slight change in the initial condition (as is characteristic for chaotic systems). As the composition map is discontinuous at points where the composition changes, it can model a sudden change in the shape of the trajectory.

---

Questions

Q1: As long as the trajectory is twice continuously differentiable and somewhat well-behaved, we do not think it can fall outside the proposed motif set. However, we agree it can fall outside the chosen composition library if, e.g., the number of inflection points is very large or infinite. Please see our response in "(A) Flexibility of motif-based semantic representation".

Q2: Response in "(B) Higher dimensions"

Q3: Response in "(C) Comparison with black-box baselines"

Q4: Response in "(E) Granularity of information"

Q5: Response in "(F) Chaotic behavior"

Q6: This paper lays the vital groundwork for direct semantic modeling with many future opportunities to extend it to even more settings, broadening its applicability. Our current focus is indeed on extending the framework by meta-motifs to accommodate some infinite compositions and to extend to multi-dimensional inputs as described in our response in "(B) Higher dimensions".

---

References

[1] Chen, R. T., Rubanova, Y., Bettencourt, J., & Duvenaud, D. K. (2018). Neural ordinary differential equations.

[2] Holt, S. I., Qian, Z., & van der Schaar, M. (2022). Neural laplace: Learning diverse classes of differential equations in the laplace domain.

[3] Lu, L., Jin, P., & Karniadakis, G. E. (2019). Deeponet: Learning nonlinear operators for identifying differential equations based on the universal approximation theorem of operators.

We have included the results in Appendix B.5, which shows the predicted trajectories and some of the property maps. The average RMSE on the test dataset is 0.019 for $S$ trajectory, 0.011 for $I$ and 0.014 for $R$. Note that the irreducible error on this dataset (caused by added Gaussian noise) is 0.01. The shown shape functions let us draw the following insights about the model:

• The time when $I$ is at its maximum $t_{\text{max}}(I)$ is on average just below $0.2$. $I_0$ has a relatively large impact on $t_{\text{max}}(I)$ by increasing it by 0.1 for very low $I_0$ or decreasing it by 0.05 for very high $I_0$. The larger the $I_0$, the faster the maximum is achieved.
• $S_0$ also has a negative impact on $t_{\text{max}}(I)$ but it is much smaller ($\pm 0.02$).
• The maximum of $I$ (denoted $I_{\text{max}}$) increases linearly with both $S_0$ and $I_0$. This time $S_0$ has slightly bigger impact ($\pm 0.1$) compared to $I_0$ ($\pm 0.04$) .
• In both $t_{\text{max}}(I)$ and $I_{\text{max}}$, the impact of $R_0$ is insignificant.
• The horizontal asymptote of $R$ increases linearly with all three initial conditions. In particular, the shape function associated with $R_0$ has a unit slope as expected.

An interesting advantage of our approach is that even though three variables describe the system, we do not need to observe all of them to fit the trajectory (similar to the pharmacokinetic example in the paper). ODE discovery methods assume that all variables are observed, which constrains their applicability in many settings.

---

(C) Comparison with black-box baselines

We have now added a comparison with three black box approaches: NeuralODE [1], NeuralLaplace [2], and DeepONet [3]. All results are now included in Table 3 and can be seen below. In nearly all settings, Semantic ODE outperforms the black-box method. This is likely because we operate in a low-data regime, where we always have fewer than 200 samples measured at most 20 times throughout the trajectory.

---

(D) Conceptual comparison with SINDy

As described above, our approach and SINDy contain assumptions on the class of solutions. We argue that often the users may have some prior knowledge about the behavior of the dynamical system (e.g., the shape of the trajectory, its asymptotic behavior, or the fact that it is always increasing at the beginning) but not necessarily about its functional form. It may also be the case that no compact closed-form ODE exists. For instance, the dynamics may depend on history (e.g., delay and integral differential equations), involve unobserved variables (e.g., pharmacological model), or empirically determined functions (e.g., stress-strain curve). Even if we are able to translate our prior knowledge into assumptions on the functional form, enforcing them may be challenging.

The library of terms in SINDy on its own is usually insufficient to determine the convexity and monotonicity of the predicted trajectories. Even a small library of terms ($x$, $t$, $1$) can give trajectories that are increasing, decreasing, convex, and concave, as well as trajectories with trend changes, depending on the exact coefficients as well as the initial condition and the time domain. So, there is a substantial difference between choosing the library of terms and the set of admissible compositions.

Dear Reviewer rHqF,

Thank you very much for such a comprehensive and insightful review. We are glad you found our work original, paper clear, and potential significance substantial. By addressing your questions and concerns, we truly believe we have strengthened our paper. We provide our responses below and point to relevant changes in the paper.

---

(A) Flexibility of motif-based semantic representation

Representation theorem

We have proved the following representation theorem that shows the generality of our used motifs.

Any function $x \in \mathcal{C}^2(t_0, +\infty)$ whose second derivative vanishes at a finite number of points (i.e., $|\{t \ : \ \ddot{x}(t) = 0\}| < \infty$) can be described as a finite composition constructed from proposed motifs.

The proof is provided in Appendix G.5.

Extension to infinite compositions

As we mentioned in the limitations, our approach is limited to trajectories with a finite composition. We have now elaborated in Appendix G.1 how this can be extended to oscillatory/periodic trajectories by meta-motifs and their associated meta-properties. The idea is to take the infinite composition (infinite sequence of motifs) and represent it using some finite representation. For instance, a trajectory $\sin(t)$ would be described as $(s_{+-c},s_{--c},s_{-+c},s_{++c})$ repeating forever, where the meta-properties may include the "frequency" and "amplitude" that may also vary with time and be itself described using compositions. We admit that this extension still assumes some underlying pattern behind the infinite sequence of motifs that can be "compressed" into a shorter representation. This is slightly related to the Kolmogorov complexity of this sequence. If the infinite sequence of motifs is algorithmically random, then motif-based representation may not be helpful. We have now made it clear in the limitations.

Dependence on a pre-defined set

Although our current implementation depends on a pre-specified set of compositions (e.g., all possible compositions of length less than $5$), ODE discovery methods also rely on the specification of the token set (e.g., a library of terms in SINDy or a list of operators and functions in symbolic regression). Both of them introduce a form of inductive bias. As we argue in Section 6.1, semantic inductive biases may be more meaningful and intuitive for users than syntactic ones. Especially when there is no prior knowledge about the functional form of the ODE that best describes the dynamical system.

---

(B) Higher dimensions

We have now elaborated on extending semantic modeling to multiple dimensions in Appendix G.2. The easiest way to model multiple trajectories (or multi-dimensional trajectories) is to model each trajectory dimension independently but conditioned on the multi-dimensional initial condition $x_0\in\mathbb{R}^M$, i.e., for $M$ trajectories, we want to fit $M$ forecasting models, $F_1, \ldots, F_M$ such that for each $m\in[M]$, $F_{m}:\mathbb{R}^M \to \mathcal{C}^2$. Each $F_m$ can still be represented as $F_{\text{traj}} \circ F_{\text{sem}}^{(m)}$, where $F_{\text{traj}}$ is the same as described in the paper and $F_{sem}^{m} : \mathbb{R}^M \to \mathcal{C}\times\mathcal{P}$. $F_{sem}^{m}$ can similarly be represented using a composition map and property maps for each predicted composition. However, these maps are no longer described by univariate functions. We currently envision the composition maps to resemble decision trees (or similar structures) and the property maps as generalized additive models (GAMs) so that the whole structure remains comprehensible. Moreover, the semantic predictors can also accommodate auxiliary variables beyond the initial conditions. For instance, the parameters of the ODEs. This approach should scale to multiple dimensions as it is very modular. Each $F_{\text{sem}}^m$ can be analyzed independently. Each property map of $F_{\text{sem}}^m$ corresponding to a particular composition can be analyzed independently. And finally, each shape function for each of the predicted properties can be analyzed independently. We believe this modularity is one of the reasons why even systems with multiple dimensions can remain understandable and, more importantly, can be edited as each change has a localized impact. In contrast, changing a single parameter in a system of ODEs may result in a change of all the trajectories in ways that may be difficult to predict without a careful analysis.

This extension deserves a separate paper, but we wanted to demonstrate a proof of concept. We implemented the property maps as we defined above (each property described as GAMs) and fitted data following an SIR epidemiological model for different initial conditions. To simplify the problem, we prespecify the composition map (we only train property maps). We assume $S$ follows $(s_{--c},s_{-+h})$, $I$ follows $(s_{++c},s_{+-c},s_{--c},s_{-+h})$, and $R$ follows $(s_{++c},s_{+-h})$.

Dear Reviewer rHqF,

We are deeply grateful for your thoughtful consideration of our responses and for increasing your score. Your feedback has been invaluable in refining our work. If there are any additional suggestions or concerns you have, we would be more than happy to address them to further improve our paper.

Thank you once again for your support and encouragement.

Kind regards,

Authors