Efficient Training of Minimal and Maximal Low-Rank Recurrent Neural Networks

We thank the reviewer for their time & valuable feedback on our paper. We especially appreciate their positive remarks on the advantages of our method compared to FORCE/BP and the clarity of our methods. Below, we address specific weaknesses and questions raised by the reviewer, and humbly request the reviewer to consider raising their score if they are satisfied with our responses. Weaknesses:

1. Performance on high dimensional systems: The reviewer is absolutely correct in highlighting this, as indeed we have not considered high-dimensional dynamical systems. However, the primary motivation for our work is to extend the theoretical foundations of the low-rank RNN literature (e.g., Mastrogiuseppe & Ostojic, 2018; Beiran et al., 2021; Dubreuil et al., 2022; Valente et al., 2022) & the “task training” RNN literature (e.g., Mante & Sussillo 2013; Sussillo & Barak 2013; Luo, Kim, et al., bioRxiv 2023), which has focused on training RNNs to perform decision-making tasks (e.g., appendix Section C of our paper). Thus, while these tasks aren’t high-dimensional (i.e typically rank 1 or 2) , they're challenging for animals to perfect & perform. Additionally, while the ODEs considered in the main paper are either 1 or 2 dimensional, in Section G.2 of the appendix, we do show application to a 3-dimensional system (the Lorenz attractor), along with a comparison to FORCE.

We would like to add, moreover, that there is no theoretical impediment to applying our framework to higher dimensional systems; the regression-based approach we have described is fully general, and thus equally applicable to systems of higher dimensionality. To illustrate this point, we will apply our method to 10, 20, and 40-dimensional systems in the revised manuscript. In fact, our marginal-likelihood based approach to fitting infinite RNNs is especially advantageous in higher dimensions, where selecting an optimal basis for a particular nonlinear function is even more critical. We feel this will substantially strengthen the paper, illustrating application to dynamical systems much higher-dimensional than those commonly considered in most neuroscience tasks.

2. Comparison to older methods (FORCE/BP): We thank the reviewer for highlighting this. However, to the best of our knowledge, the majority of the literature we build on (referenced in point above) uses these frameworks for training low-rank RNNs, which is why we focus our comparison on them.

Questions:

1. v=1, poor approximation and \hat{n} computation: We apologize for the confusion. For v=1, our approximation is nearly perfect (please see Fig 1 panel E,F). The approximation is poor in the case of no inputs (v=0), as highlighted in Fig 1 Panel C,D. Lastly, if Iv is obtained from a Gaussian distribution, the computation remains linear in ‘n’, thus maintaining the same Least Squares solution presented in equation 8.

2. Approximation of g(z) in high-dimensional settings: Thank you for pointing this out. The approximation is achieved in the same way as the examples presented in the paper. We briefly discuss how this approximation is achieved in the following two cases:

   2.1. If the complete ODE is known: In the case of a high-dimensional dynamic, if the complete ODE is known, g(z)’s approximation would simply be the evaluation of the different dimensions of ‘g’. Thus each component of g(z) would be stacked to obtain the complete output. A complete description of this can be seen through our exposition of the 2-d limit cycle in Sec 3.1 (additionally Sec G.2 of the appendix extends this to the 3-d case as well).

   2.2. Only trajectories are available (target matching scenario): Alternatively, in the case of target matching, where trajectories from the ODE (g(z)) are available (eg see Sec 4 of the paper); we can approximate g(z) via finite differencing (Euler method). Thus $g(z) \approx \frac{z_t - z_{t-1}}{dt}$, where $dt$ represents the Euler bin.

Thus the framework presented is applicable to any dimensional ODE, however the corresponding size of the target and design matrices would increase.

3. Application of method to tasks and broader applicability: Thank you for this question. The goal of our framework isn’t restricted to analyzing time-series data; rather it is designed to train a low-rank RNN to replicate a target dynamical system. Keeping this view in mind, our framework goes beyond simply approximating g(z), but rather provides an interpretable and efficient way to train low-rank RNNs.

   3.1. Specifically as highlighted in the literature (Elia,et al.,2021; Barak 2017; Vyas, Saurabh, et al. 2020), this trained model becomes a proxy for the real system. This allows us to do experiments such as identifying minimal neurons needed, neuron ablations, perturbation analysis etc which might not be possible to carry out in the true biological or physical system (or might be costly/time intensive). Thus, such a proxy would accelerate scientific understanding of such systems. However, such models have typically been seen as black boxes, and we believe our work provides stronger theoretical and practical insights into the internal representations (and learning dynamics) of these networks.

   3.2. Second, we also present an active learning strategy which provides a principled way to identify informative data points, thereby reducing experimental costs.

We thank the reviewer for their time and positive assessment on the originality and clarity of our work. We’re extremely grateful that they recognize the value of our contributions and share our excitement about the future directions this line of research can take! Below we address the specific weaknesses and questions raised by the reviewer:

Limited applicability to neuroscience: We appreciate the reviewer for giving us a chance to address this, and apologize for our lack of discussion along these lines in the paper. We believe our work provides an alternative way to accelerate neuroscientific research. Specifically, typically such analyses of task dynamics is done in a data-driven manner. That is, task-trained RNNs interpret governing dynamics of specific tasks linked to brain computation (e.g., Kanitscheider & Fiete 2017; Turner, E, et al. 2021; Luo, Kim, et al. 2023). While this approach has numerous advantages, our lack of understanding of how such networks embed solutions leads to degeneracies in the network/task solution space therefore often leading to incorrect/evolving conclusions. It is for this reason, we propose our work as an alternative — “hypothesis-driven” analysis. To elaborate, we propose designing ODEs to solve tasks (or have access to trajectories from such ODEs), which our method then embeds into the low-rank RNN. As highlighted in our work, this reverse first process allows us to make interesting and concrete explanations regarding the network connectivity, the minimal number of neurons, ideal network initializations and a first step to closed loop experiments via active learning. We do however agree with the reviewer, that this is limited to tasks for which experimentalists/computational collaborators can effectively hypothesize such underlying low-dimensional dynamics. We will add a discussion on this in the final version of the paper.

Below we address specific weaknesses highlighted by the reviewer:

1. Fairly idealized, low-dimensional and noise-free examples: The reviewer is correct in their intuition, as the latent dynamic dimension increases the required basis dimensionality also increases. While the examples in the paper are mostly 1 or 2D; we have included the 3D Lorenz attractor example and comparison of our RLS method with FORCE in the appendix (Section G.2.1). Note, in this case 16 neurons (for our RLS method) suffices for almost perfect reconstruction. To the reviewers second point on noise-free examples, we have included new analyses on a noisy target from a sum of four sinusoids (as done in the original FORCE paper). As can be seen from our analysis, the reviewer’s intuition is correct, i.e our RLS method converges faster and with far fewer neurons ( details under response to Q2 below).
2. Fitting low-rank RNNs to known ODEs/ PCA on him dim noisy data: We thank the reviewer for highlighting this. On fitting low-rank RNNs to known ODEs – we point the reviewer to our response to the main weakness highlighted by them (Limited applicability to neuroscience). Once again, while we agree designing such an ODE for a potential task might not always be possible, we believe that in the cases that it is (such as typical 2AFC decision making/ cognitive tasks studied in this literature), our framework provides a rigorous and insightful solution.
   • Point on PCA: we agree with the reviewer, our goal was simply to illustrate that our framework can also fit such trajectories if desired. However, we believe embedding the latent dynamics into a trained low-rank RNN offers capabilities beyond what PCA (since it’s a linear dimensionality reduction) alone provides such as understanding the influence of perturbations etc.
3. Minimal size RNN: The reviewer is correct in their assessment, the main advantage of OMP/COMP is from a computational perspective of fitting small networks.
4. Infinite low-rank RNNs: We appreciate the reviewer for highlighting this! To address this comment we’ve conducted new analysis, please see our result details under Q5 below.
5. Fit to real data: The reviewer’s assessment is correct. Our goal was to present interpretability from a computational perspective. I.e our method provides a geometric/algebraic description on how each neuron (in the artificial network) is a non-linear basis function of the ODE of interest. We agree there might be intriguing neuroscientific interpretations of this framework, however exploring these directions is left for future work.

Questions:

1. Training time RLS, BPTT (only n): We simulated a target signal (T=1000s, dt =0.1s, thus 10,000 samples) of a sum of four sinusoids. We then added independent gaussian noise $\mathcal N(0,0.05)$ at each time step of the trajectory, thus obtaining a noisy target as seen in FORCE. While the output is a 1-d signal, we used the trajectory and two cumulative sums as our target (thus a rank 3 network). We trained low-rank networks with BPTT only allowing ‘n’s weights to be updated. Networks were trained between 5-8 epochs (i.e to achieve similar error as our method) and average running time is reported below across 5 seeds. Note, our RLS method converges faster.
   | Number of Neurons| Method | Time (s) | |----------|----------|----------| | 8 | RLS (rank 3) | 0.07 | | 8 | BPTT ( rank 3) | 105.28 | | 8 | BPTT ( rank 4) | 83.29 | | 16 | RLS (rank 3) | 0.12 | | 16 | BPTT ( rank 3) | 70.01 | | 16 | BPTT ( rank 4) | 59.35 | | 64 | RLS (rank 3) | 0.28 | | 64 | BPTT ( rank 3) | 138.82 | | 64 | BPTT ( rank 4) | 128.69 | | 128 | RLS (rank 3) | 0.64 | | 128 | BPTT ( rank 3) | 153.50 | | 128 | BPTT ( rank 4) | 28.92 | | 256 | RLS (rank 3) | 1.41 | | 256 | BPTT ( rank 3) | 131.4 | | 256 | BPTT ( rank 4) | 140.03 | | 512 | RLS (rank 3) | 8.34 | | 512 | BPTT ( rank 3) | 130.0 | | 512 | BPTT ( rank 4) | 130.93 | | 1024 | RLS (rank 3) | 41.21 | | 1024 | BPTT ( rank 3) | 129.89 | | 1024 | BPTT ( rank 4) | 128.98 |
2. RLS v FORCE v BPTT: We use the same noisy sum of four sine-waves target discussed above. RNNs are trained with BPTT and all weights are updated. To ensure a consistent comparison, we report the training MSE (averaged across 5 seeds) after one epoch (a single pass over the 10,000-step sequence). Both RLS and FORCE achieve near-perfect performance in a single pass (note our method has highest performance), while BPTT typically requires 5–8 epochs to reach comparable errors. Here, we report BPTT results after one epoch to highlight the difference in time efficiency (as we can’t include images of training loss curves). | Number of Neurons| Method | Train Error (MSE) | Std Error | |----------|----------|----------|----------| | 8 | RLS (rank 3) | 0.012 |0.006 | | 8 | FORCE | 0.525 | 0.064 | | 8 | BPTT (rank 3) | 0.107 | 0.106 | | 8 | BPTT (rank 4) | 0.087 | 0.0919 | | 16 | RLS (rank 3) | 0.012 | 0.006 | | 16 | FORCE | 0.530 | 0.029 | | 16 | BPTT (rank 3) | 0.047 | 0.048 | | 16 | BPTT (rank 4) | 0.070 | 0.0768 | | 64 | RLS (rank 3) | 0.011 | 0.006 | | 64 | FORCE | 0.275 | 0.056 | | 64 | BPTT (rank 3) | 0.016 | 0.003 | | 64 | BPTT (rank 4) | 0.017 | 0.006 | | 128 | RLS (rank 3) | 0.011 | 0.006 | | 128 | FORCE | 0.127 | 0.058 | | 128 | BPTT (rank 3) | 0.016 | 0.003 | | 128 | BPTT (rank 4) | 0.016 | 0.003 | | 256 | RLS (rank 3) | 0.011 | 0.006 | | 256 | FORCE | 0.024 | 0.006 | | 256 | BPTT (rank 3) | 0.023 | 0.011 | | 256 | BPTT (rank 4) | 0.019 | 0.009 | | 512 | RLS (rank 3) | 0.011 | 0.006 | | 512 | FORCE | 0.014 | 0.002 | | 512 | BPTT (rank 3) | 0.022 | 0.013 | | 512 | BPTT (rank 4) | 0.020 | 0.009 | | 1024 | RLS (rank 3) | 0.011 | 0.006 | | 1024 | FORCE | 0.012 | 0.001 | | 1024 | BPTT (rank 3) | 0.019 | 0.004 | | 1024 | BPTT (rank 4) | 0.017 | 0.003 |
3. Qualitative recovery: Yes, we achieve MSE within 10^-3, and the flow-fields when reproduced do also qualitatively match.
4. Sparsity of n: This is an interesting performance comparison which we will work on including in our final version of the paper! Thank you for bringing this to our attention.
5. Optimized basis: As per the reviewer’s suggestion we have also included new analysis on finite sized low-rank RNNs trained with standard normal vs optimized basis for the ODE depicted in Fig 5 of our paper. We re-frame this problem as a target matching problem. Specifically, we used a set of four starting locations given by $[7.92, 4.72, 2.3, -7.36]$. We then rolled out the dynamics (as per the ODE in Fig 5) using an Euler integration for a total of 1000 time-steps (with dt=0.1). Each of these trajectories converged to a specific stable fixed point. We then compared low-rank RNNs (rank=1) trained via our RLS method (trained for 1 epoch) and those trained with BPTT (trained for 5 epochs) on standard normal initialized basis and the optimal basis distributions found via our GP framework. Below we summarise the Train MSE. Note, in both finite size RNNs trained with our RLS and those trained with BPTT, the optimized basis results in lower train errors. | Number of Neurons| Method | Train MSE | |----------|----------|----------| | 2 | RLS | 0.847| | 2 | RLS (opt) | 0.125 | | 2 | BPTT | 6.519 | | 2 | BPTT (opt) | 0.211| | 5 | RLS | 0.021 | | 5 | RLS (opt) | 0.001 | | 5 | BPTT | 0.508 | | 5 | BPTT (opt) | 0.199| | 10 | RLS | 0.041 | | 10 | RLS (opt) |0.001 | | 10 | BPTT | 0.250 | | 10 | BPTT (opt) |0.187 | | 50 | RLS | 1e-3| | 50 | RLS (opt) | 1e-3 | | 50 | BPTT | 0.239 | | 50 | BPTT (opt) | 0.1907| | 100 | RLS |1e-4 | | 100 | RLS (opt) | 1e-5 | | 100 | BPTT | 0.200 | | 100 | BPTT (opt) | 0.188|

Additional Questions: We thank the reviewer for bringing these to our attention. These suggestions seem very interesting! While not in the scope of this rebuttal, we will include these analysis in future work!

We thank the reviewer for their time and positive assessment of our work. We especially appreciate their constructive suggestions on improving our work.

Below we address specific questions raised by the reviewer:

1. Optimal error rates and number of neurons: We’re grateful for being pointed to this literature, it certainly relates to our empirical results on OMP/COMP, and active learning. We agree with the reviewer that connecting our results with the theoretical literature on optimal error rates will certainly strengthen the paper, and we will add a discussion (and citations) to this in our final version!

2. Clever kernel choices: Yes certainly! Thanks for pointing this out, we will include a reference to this.

3. Effect of sampling (z,I): Our active learning framework builds on previous results in bayesian optimal/adaptive experimental design. Specifically, we cast the problem in a linear-gaussian setting and use the Bayes A-optimal criteria (i.e maximizer of predictive posterior variance) to choose the next point. Within this setting, we can compute the minimum number of samples required to obtain a desired MSE; this bound depends on the dimensionality of the kernel function over the support of our system (i.e., the eigenspectrum of $\Phi(x)^T \Phi(x)$ for a grid of ‘x’ values over this support) and the level of (Gaussian) observation noise in each sample (see Chaloner & Verdinelli 1995). We will include an explicit discussion on this in our final version!

   To further validate our empirical results we also performed additional analysis comparing our active learning method to random sampling through for the ODE presented in Fig 6. Note, our method achieves significantly lower error rates much quicker. Below we summarize our results (the first data point is selected randomly for our method as well for fair comparisons), MSE is computed on average over $5$ random seeds: | Number of data-points | Sampling Method | MSE (True ODE vs Fitted) | |----------|----------|----------| | 1 | Random | 0.814 | | 1 | Ours | 0.814| | 2 | Random | 0.820 | | 2 | Ours | 0.806 | | 3 | Random | 0.478 | | 3 | Ours | 0.010 | | 4 | Random | 0.012 | | 4 | Ours | 0.006 | | 5 | Random | 0.011 | | 5 | Ours | 0.006 | | 6 | Random | 0.008 | | 6 | Ours | 0.005 | | 7 | Random | 0.007 | | 7 | Ours | 0.003 | | 8 | Random | 0.005 | | 8 | Ours | 0.002 | | 9 | Random | 0.002 | | 9 | Ours | 0.0009 | | 10 | Random | 0.002 | | 10 | Ours | 0.0008 |

4. Initial parameter initializations and minimal rank bias: We deeply apologize for the lack of details, which we will include in the final version of the paper. For the initial comparisons in Fig 3 (for our method and BPTT), we drew $m, I$ from $\mathcal N(0,1)$. FORCE needed initializations $\mathcal N(0,\frac{1}{d})$, where $d =$ number of neurons. For the original set of experiments, the results presented averaged over $3$ seeds of initialization.

   • Minimal Rank Bias: We thank the reviewer for bringing this to our attention. We’re currently working on training full-rank networks to reproduce a noisy target ((which requires a minimum of rank 3 per our framework). We will check how different initializations correlate with the effective rank of the recurrent matrix (similar to Bordelon et al. (2025)). We'll be sure to include this comparison in the final version of the paper!

5. Discussion on embedding known ODEs: We apologize for our lack of discussion on this topic and will include a section in our camera-ready version. Here we provide a brief explanation on this. We believe our work provides an alternative way to accelerate neuroscientific research. Specifically, typically such analyses of task dynamics is done in a data-driven manner. That is, task-trained RNNs interpret governing dynamics of specific tasks linked to brain computation (e.g., Kanitscheider & Fiete 2017; Turner, E, et al. 2021; Luo, Kim, et al. 2023). While this approach has numerous advantages, our lack of understanding of how such networks embed solutions leads to degeneracies in the network/task solution space therefore often leading to incorrect/evolving conclusions. It is for this reason, we propose our work as an alternative — “hypothesis-driven” analysis. To elaborate, we propose designing ODEs to solve tasks (or have access to trajectories from such ODEs), which our method then embeds into the low-rank RNN. As highlighted in our work, this reverse first process allows us to make interesting and concrete explanations regarding the network connectivity, the minimal number of neurons, ideal network initializations and a first step to closed loop experiments via active learning.

We thank the reviewer for their time and enthusiastic assessment on the clarity, originality, and significance of our work. We’re pleased that they found our paper clear and well-organized and appreciated both our new theoretical and practical insights. We’re especially encouraged by their recognition of our effort to bridge multiple frameworks in RNN learning theory and to provide new insights supported by sound experiments!

Below we address specific questions raised:

1. Difference from NEF: We apologize for the lack of clarity on this front. Our presentation is exactly that of NEF in the case of the entire ODE being known. We would like to highlight, although the NEF presentation is indeed low-rank, we haven’t seen any explicit literature highlighting this (to the best of our knowledge). Additionally, we believe our online RLS approach is the first formulation to adapt NEF for the target matching case. While this appears straightforward, we believe that we’re the first to present this. Similarly, our other novel contributions such as finding the smallest RNN (OMP/COMP), linking the framework to the GP formulation and active learning are all novel connections to the NEF methodology. We will include a detailed discussion on this in our final version of the paper.
2. lines 99-100: We meant to highlight a case where a task/target may need more latent dimensions than the output signal. For instance, a 1-d periodic signal would need at least 2-d latents to represent oscillations. We apologize for the confusion, and will simplify our text to more clearly indicate this.
3. lines 128-129: This refers to our discussion on networks with the relu non-linearity. Here, we highlight, similar to tanh networks in the absence of per neuron inputs (or biases), the basis functions are piecewise functions about origin only, and can thus only universally express functions that are not offsetted from the origin (Fig 7 in the appendix).
4. Equations 8, 10: We thank the reviewer for bringing this to our attention, we’ll clean up our notation in the final version of the paper!
5. Variance of m,I: We deeply apologize for the lack of details, which we will include in the final version of the paper. For the initial comparisons in Fig 3 (for our method and BPTT), we drew $m, I$ from $\mathcal{N}(0, 1)$ and for FORCE we used $\mathcal{N}(0, \frac{1}{d})$, where $d=$ number of neurons.
6. Tasks presented are too simple: We agree with the reviewer that the tasks presented in the paper are fairly simple. Our goal was to present simpler toy tasks in the main paper to highlight our theoretical advancements and to develop intuition for our presented results. Below we discuss various additional/more complex tasks that we have included.
   • 3-d example, contextual tasks: The reviewer is absolutely correct in highlighting this, as indeed we haven’t considered high-dimensional dynamical systems. However, the primary motivation for our work is to extend the theoretical foundations of the low-rank RNN literature (e.g., Mastrogiuseppe & Ostojic, 2018; Beiran et al., 2021; Dubreuil et al., 2022; Valente et al., 2022) & the “task training” RNN literature (e.g., Mante & Sussillo 2013; Sussillo & Barak 2013; Luo, Kim, et al., bioRxiv 2023), which has focused on training RNNs to perform decision-making tasks (e.g., appendix Section C of our paper). Thus, while these tasks aren’t high-dimensional (i.e typically rank 1 or 2) , they're challenging for animals to perfect & perform. Additionally, while the ODEs considered in the main paper are either 1 or 2 dimensional, in Section G.2 of the appendix, we do show application to a 3-dimensional system (the Lorenz attractor), along with a comparison to FORCE. We would like to add, moreover, that there is no theoretical impediment to applying our framework to higher dimensional systems; the regression-based approach we have described is fully general, and thus equally applicable to systems of higher dimensionality. To illustrate this point, we will apply our method to 10, 20, and 40-dimensional systems in the revised manuscript. In fact, our marginal-likelihood based approach to fitting infinite RNNs is especially advantageous in higher dimensions, where selecting an optimal basis for a particular nonlinear function is even more critical. We feel this will substantially strengthen the paper, illustrating application to dynamical systems much higher-dimensional than those commonly considered in most neuroscience tasks.
   • Noisy/Complex Target: During this rebuttal period we have also conducted new analyses on a more complex: noisy sum of four sine-waves target signal. Specifically: we simulated a target signal (T=1000s, dt =0.1s, thus 10,000 samples) of a sum of four sinusoids. We then added independent gaussian noise (N(0,0.05)) at each time step of the trajectory, thus obtaining a noisy target as seen in FORCE. While the output is a 1-d signal, we used the trajectory and two cumulative sums as our target (thus requiring a rank 3 network), for our RLS method. We compared our method to low-rank RNNs trained with BPTT (where all weights were updated), as well as FORCE networks (note: all weights were sampled from $N(0,\frac{1}{\sqrt d})$ as referenced by FORCE). To ensure a consistent comparison, we report the training MSE (averaged across 5 seeds) after one epoch (a single pass over the 10,000-step sequence). Both RLS and FORCE achieve near-perfect performance (RLS has best performance) in a single pass, while BPTT typically requires 5–8 epochs to reach comparable errors. Here, we report BPTT results after one epoch to highlight the difference in sample efficiency (as we can’t include images of the complete training loss curve).

• Comparison with optimized v un-optimized basis: Finally, we have also included new analysis on RNNs trained with standard normal vs optimized basis. We reframe this problem as a target matching problem. Specifically, we used a set of four starting locations given by $7.92, 4.72, 2.3, -7.36$. We then rolled out the dynamics (as per the ODE in Fig 5) using an Euler integration for a total of 1000 time-steps (with dt=0.1). Each trajectory converged to a specific stable fixed point. We then compared low-rank RNNs (rank=1) trained via our RLS method and those trained with BPTT (for 5 epochs) on standard normal initialized basis and the optimal basis distributions found via our GP framework. Below we summarise the Train MSE. Note, in both finite size RNNs trained with our RLS method and those trained with BPTT, the optimized basis results in lower train errors.

7. Optimization of cov, ref [3], limitations of a chosen basis: We appreciate the reviewer for giving us a chance to address this. The choice of $m,I$ matches the distribution which we originally used in our comparison. As per our understanding, ref [3] and follow up work (Dubreuil, et al. 2022; Valente, A. et al. 2022 ) suggests a joint normal distribution with $m,I$ (not necessarily with $0$ mean as we have used in our implementations and theoretical results).
   • Mean $0$ distributions: While for this version we have worked with $0$ mean distributions (as this cleanly relates to the arcsin kernel), including the mean as a kernel hyper-parameter has been left to future work. We will clarify our claim in the paper, thanks for pointing this out!
   • Limitations of a choses basis: The reviewer brings up a very interesting point on when correlations with $m, I$ are needed. While we agree this is an intriguing question (one which we do not answer in this work), we provide an empirical framework to optimize such a basis, to be as effective as possible in approximating an ODE. We do however agree with the reviewer and will include an explicit discussion on this in the final version of our paper!

Lastly, we believe that our framework offers the same capabilities (with smaller, faster trained and more interpretable networks) as that of FORCE/NEF or BP (in the case of target matching tasks), as demonstrated in our results in Fig 2, appendix G.2 and our new analysis on a noisy target trajectory.