{'fraction_non_alphanumeric': 0.03227973736448313, 'fraction_numerical': 0.011146739960299282, 'mean_word_length': 4.468603874415497, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 2, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Inferring the dynamics of oscillatory systems using recurrent neural networks

We investigate the predictive power of recurrent neural networks for oscillatory systems not only on the attractor, but in its vicinity as well. For this we consider systems perturbed by an external force. This allows us to not merely predict the time evolution of the system, but also study its dynamical properties, such as bifurcations, dynamical response curves, characteristic exponents etc. It is shown that they can be effectively estimated even in some regions of the state space where no input data were given. We consider several different oscillatory examples, including self-sustained, excitatory, time-delay and chaotic systems. Furthermore, with a statistical analysis we assess the amount of training data required for effective inference for two common recurrent neural network cells, the long short-term memory and the gated recurrent unit.

Inference of most dynamical properties of any system is typically best done with an active experiment, meaning that one has the power to repeatedly manipulate the system state in controlled conditions in order to isolate the desired measure, e.g. car technical inspection, where the examiners manipulate the car into specific situations in order to asses its safety on the road. Often however, one only has access to passive observations. This can be due to a number of reasons, for example, the system can be very big, like when studying the dynamics of our planet, or the system can be delicate, like when studying the phisiology of the human body. The algorithms that attempt to distil dynamical measures from passive observations commonly make assumption on how the observations were collected and typically require very long observations. Imagine, for example, assessing the car's capabilities only from observing routine trips to work. Here we propose a conceptually simple scheme, relying on the now well established artificial neural networks. In particular, we use recurrent neural networks which have established themselves in timeseries forecasting, text generation etc. and train the networks to mimic the system dynamics, allowing us to then perform an active experiment on the trained model. We test this on several oscillatory systems and measure their characteristic properties, such as, bifurcations, dynamical response curves, characteristic exponents and compare them to the measures form the original system. a) Electronic mail: rokcestn@uni-potsdam.de, r.cestnik@vu.nl

I. INTRODUCTION
Oscillatory systems can be found in all fields of natural science: in optics 1 , electronics 2 , chemistry 3 , biology 4 , climatology, life science, etc. Oscillations are present at all scales, from biological cells like neurons 5 , to organs like the heart 6 , to oscillations spanning the entire organism such as the circadian rhythm 7 and the menstrual cycle to many years, as for some planetary motion.
For technical systems, the classical modeling of oscillatory systems consists of hard thinking which terms are involved, followed by tests of the model concluded. For low-dimensional systems, this works very well; however the aforementioned examples are all complex, high-dimensional and coupled to their surrounding, like the brain which consists of many coupled neurons 8 , climate models 9 or fluids, whose dynamics has been a major driving force for the investigation of periodic motion, synchronization of oscillatory systems 10 , period doubling bifurcations, and chaotic oscillations. A detailed modeling of such systems is very hard, and with increasing computer power existing methods to infer dynamical systems from measurements are easier to realize 11-13 .
For high-dimensional systems, one has to either measure with many channels or apply embedding methods 14,15 , or combine both ways to analyze the data. If a system is truly periodic, the system lives on a one-dimensional manifold and the resulting system of equations may be modeled by a simple two-dimensional, possibly nonlinear system. If weakly perturbed, under certain assumptions 16 , the system moves close to the original orbit. Such perturbations may originate from another oscillator, a network of oscillators or elsewhere from the environment. However, if the system is close to a bifurcation and the corresponding parameter is perturbed externally or parametrically, the system may undergo dramatic changes in its dynamic. A bifurcation is, however hard to predict for heuristic models, whereas it is easy if equations are known. Under this point of view previous approaches using symbolic regression methods 12,17,18 proved successful. Heuristic methods like liquid state machines, echo state networks, or various types of artificial neural networks 19-21 perform very well in predicting the dynamics. However, few studies are known for particular aspects of oscillatory systems inferred from time series. Here, we investigate several oscillatory models under perturbation, as they may occur in real measurements. Our focus is on the inference of bifurcation behavior and possibly chaos, even if not all of the parameter variation is included in the measurement.
Technically, in all of the above methods one typically assumes a priori a model or a class of models, sets an optimization criterion (e.g. least squares) and optimizes model parameters or its functional constituents for nonparametric methods. Mathematical aspects are most often left aside, e.g. basic assumptions on the existence of solutions, robustness under perturbation, in particular for heuristic methods. Here we utilize the widely used artificial neural networks (ANNs). An ANN has several hyperparameters such as the actual topology of the network, the activation function, or the learning rate, and is in general very pliable toward many different tasks. Since we consider time series, we investigate the capacities of recurrent neural networks 22 (RNN). Due to loops in their connectivity they retain past information, i.e. they inherently possess a memory, similar to embedding. Further, they have directed connection which seems to be particularly successful in speech recognition 23 , text generation 24 and machine translation 25 , where a forward-oriented semantic is present. The aim of this study is to evaluate how suited RNN are for model oscillatory systems under the aspect of parameter change and perturbances. In this way, the inferred model of the oscillator can be probed via changing the perturbative signal, effectively allowing the performance of an active experiment.
The article is structured as follows: in section I A we refer to relevant related works and briefly recall the RNN method. in section II A we introduce the dynamical inference setup and training scheme. We then present numerical tests with example systems in section II B where we compare signal reconstruction and other observables such as phase response curve 26 (PRC) and the maximum Lyapunov exponent 27 . The different example systems are chosen as representatives of different mechanisms giving rise to oscillatory behavior; specifically, self-sustained and excitatory oscillations, time-delay induced oscillations and chaos. We continue by presenting numerical tests on data requirement for successful inference in section II C where we compare the inference quality for different lengths of time-series used for training. We present the methods used in more detail in section III and finally, discuss the novelty, limitations and generalizations of our approach in the discussion section IV.

A. Previous work
In this paragraph we chronologically go through works related to this paper. In Ref. 28 the author uses a RNN for learning state space trajectories. In Ref. 29 the authors show that any trajectory generated by a finite-dimensional dynamical system can be effectively represented as with neural network. In Ref. 30 the authors model a dynamical system with a perturbation using a RNN. In Ref. 31 the authors use feed-forward neural networks 32 (FNN) to model dynamical systems. They feed in delayed values of one variable as well as a control parameter as inputs and train the network for one step predictions. The approach works well and they reproduce bifurcation diagrams of several example dynamical systems. In our approach we train RNNs for one step prediction where the input consists of several time-delayed values of one or more variables as well as an arbitrary number of perturbative signals. The past values of one or more variables contain information of the topology of the attractor of the complete system according to the Takens' delay embedding theorem 15 . The RNN topology prioritizes more recent values over older for the next prediction therefore we believe it more suitable for time-series prediction and demonstrate its efficiently throughout this paper.

B. Recurrent neural networks
Artificial neural networks is nowadays a relatively broad term as many different network topology classes are commonly used for different types of problems. The simplest class of ANNs are the feed-forward networks, they have directed connections between subsequent layers without any loops, effectively allowing the information to flow in only one direction -forward. The slightly more general class are the recurrent neural networks, they can have loops in their connectivity, which can result in internal state memory. Different RNNs then differ in the fine architecture of the basic cells, the order and type of logical operations. In this work we apply two commonly used cells, the long short-term memory cell 33 (LSTM) and the gated recurrent unit 34 (GRU). LSTM was constructed first in an attempt to deal with long term dependencies and GRU emerged as its faster simplification. The software implementation was accomplished with the help of TensorFlow 35 and Keras 36 .

A. Inference scheme
Consider a general dynamical system˙ x(t) = f ( x) ∈ R N , perturbed by an external perturbation p(t) ∈ R N . Suppose we have measured the timeseries of n x ≥ 1 state variables x = (x 1 , x 2 , . . . , x nx ) as well as the n p timeseries of the perturbation p = (p 1 , p 2 , p 3 , ...). The question we investigate now is if it is possible to recover both the autonomous dynamics of the system x and the system's response to the perturbation using RNN. To answer this question, we first clarify assumptions we use. Firstly, we assume the perturbation to be additive, small in amplitude and high in frequency, i.e. it does not change the general dynamical state of the system on average over a period which is significantly smaller than the smallest period of the unperturbed system. Secondly, the attractor of the unperturbed system is supposed to live on an M-dimensional manifold. Embedding theorems ensure that we need 2M + 1 variables to reconstruct the system, but less may do as well. Without perturbation we can only recover the dynamics on the attractor, but with a perturbation the phase space around the attractor is explored and we have a means to infer the neighboring phase space, too.
As implementation of the RNN we put in historical values of x(t), p(t) and return the time-evolved state x(t+∆ t ). The above explained "unrolling" of the RNN is equivalent to a corresponding time-delay embedding, where the embedding delay is constant. In RNN language, one says the network is "unrolled" by a finite number of rolls. As stated before, according to Takens theorem, the appropriate number of rolls R ≥ 2M +1 where M is the dimension of the attractor. The sampling ∆ t must be smaller than the smallest time scale which occurs in the system (or which we may want to include in our modeling). Heuristically, one can say that the time-resolution ∆ t should be chosen fine enough to see the details of interest.
Given the resolution and number of rolls we can begin to "train" our model, i.e. to start a loop for the statistical inference method: At each training step the network "learns" the possible relation: We use a least-squares optimization criterion xe− x x− x (where · stands for the mean) to determine qunatitatively how well the estimates x e match the true values x. Hereafter, we use estimated and modeled as synonymous.

B. Examples
We test our scheme on three representative model systems: a time-delay oscillator, an excitatory system and one chaotic oscillator. The validation test consists of comparing the modeled signal with the original when presented with data never seen in training. As an important measure for oscillations we estimate the phase response curve 26 (PRC) and the maximal Lyapunov exponent 27 for comparing the predictive power the model has in a dynamical systems context.
For all examples shown in this paper we use a network with 1 hidden layer of 32 nodes and 36 rolls. We use tanh activation for all but the output layer, where we use linear activation so that a continuous signal can be produced. There are two common cells used in RNN: the long shortterm memory cell 33 (LSTM) and the gated recurrent unit 34 (GRU). We tested both for the systems in this study; as a result we found that GRU performed poorly, hence all results shown are for LSTM models, cf. Sec. II C for a comparison of the two cells. To generate the data we first simulate the perturbation signal using the stochastic Euler integrations scheme, and then integrate the dynamical equations with fourth order Runge Kutta. We use a sufficiently small timestep and then re-sample the signals to an appropriately lower time reosulution to create the network training input. The resolution is chosen such that 36 points (the number of considered historical values R) corresponds to 1 natural period of the oscillator. In the case of chaotic oscillators, this was computed as the average period, in the case of excitatory systems the time needed to return from the excited state to the fixed point was used, for the timedelay oscillator we can use the known, exact period.

Roessler oscillator -phase response curve, bifurcation diagram and Lyapunov exponents
For our first test we use the Roessler system 37 , because it exhibits many different regimes from periodic oscillations with one or more frequencies, together with chaotic behavior by varying just one parameter b, cf. Fig. 2a for the bifurcation diagram. The corresponding equations, including the perturbation p read: with parameters a = 0.2 and c = 5.7. To explore phase space, we can vary b through a constant term in the perturbation p(t). For the first test we set b = 2.0, such that the system shows one frequency. In the following we use a stochastic perturbance, to generate different noise distributions we use an Ornstein-Uhlenbeck 38 process as basis: where ξ is Gaussian white noise ξ(t)ξ(t ) = δ(t − t ), and = 0.5 and τ = 5.0 are the amplitude and correlation time of p: q(t)q(t ) = 2 e −(t−t )/τ . Now, we set p(t) = q(t) and feed both the signal x(t) and the perturbation signal into the network as described in section II A. The timeseries length corresponds to 1000 natural periods, which is presented to the network during 500 training epoch in batches of 100 time points with resolution ∆ t = 0.17 (the time step used for the integration is significantly smaller). The network is trained using stochastic gradient descent 39 with learning rate 0.005.
The network learns to reproduce the dynamics to a mean deviation of 2.5×10 −2 (for the time window in Fig. 1), such that the reproduced signal is visually indistinguishable from the one generated with Eqs. (2). This holds true for both the perturbed signal -where p(t) is fed to both, the network and the equations, as well as for the unperturbed signalwhere the network is fed p(t) = 0.
Can we use the inferred network for more than just mimicking a signal, e.g. to study dynamical regimes? We want to study this scenario in probing the network for dynamical responses to stimuli. Since the system in question is a self-sustained oscillator it is natural to estimate its PRC, cf. Sec. III B. The comparison of the estimate obtained from the RNN and the the true one, is displayed in see Fig. 1. The coincidence is very good, up to mean deviation of 0.1 in the entire phase range [0, 2π). Indeed, this can be an effective method of inferring the PRC from data, cf. [40][41][42] . 1. a: The training signal generated with Eqs. (2) in red and the RNN reproduced signal in green. Both have the same input realization, depicted with gray (scaled for being visually comparable to the signal). The vertical black line marks the beginning of the forecast. b: The true phase response curve of system (2) in red and the one inferred from the RNN in green.
As stated above, our aim is to study in how far a RNN reconstruction of a system can reflect different dynamical regimes. We perform this study with the Roessler system and one trained RNN. The bifurcation parameter in the Roessler system (2) is the constant b. since we apply external driving in the y-component, we can incorporate this variation in p(t) as a constant part. This way our external driving is as well a parametrical driving of the system through the parameter b. One way to scan a whole range of parameters effectively through p is a transformation to positive values such that where q(t) is the input described in Eq. (3). It yields a process with a log normal distribution: P (p) ∼ 1 p exp − 2 log(p) + log(2) 2 . Such input spans a wide range of values, effectively introducing different regimes of our system to the network, see Fig. 2a for the input probability distribution with respect to b bifurcation (grey shaded region in the background). If we let b = 0 we can study the trained RNN with respect to this kind of input. The idea is that then, that it effectively learns to mimick the regimes corresponding to different values of b, which we can invoke via the offset of the input p(t). For this study we use a longer timeseries corresponding to 10000 natural periods, we train the network over 1000 epochs.
As a result, we find that the network reproduces the signal perturbed by a realization of Eq. (4) well. Furthermore, we can estimate the bifurcation diagram from the network by feeding it different values of constant input and observing the stationary signals, see Fig. 2a. In the value range of the input (4) the diagram obtained from the RNN matches the true one closely. It reproduces simple oscillatory regimes, chaotic regimes and the period doubling bifurcation. Throughout the range of b the natural frequency (average frequency in case of chaos) matches the true one closely, with mean deviation of 5 × 10 −2 .
In chaotic regimes the maximum Lyapunov exponent 27 is an important measure as it quantifies the divergence of nearby trajectories in time. To understand the predictive power of the RNN, we confront true value with estimate in Fig. 2b. This graph is accomplished by long time observation of the evolution of two nearby states, while re-scaling their difference to prevent them from diverging far from each other, see section III C for further details.  (2) in red and the one inferred from the RNN in green. The probability density of the input is underlaid in gray (scaled for being visually comparable to the diagram). b: The Lyapunov exponents of system (2) in red and the ones inferred from the RNN in green. The true bifurcation diagram is underlaid in gray (scaled for being visually comparable). Note that the bifurcation diagram as well as the Lyapunov exponent range were reproduced with a single RNN trained on correlated noise input, Eq. (4).
If we train a RNN with the full range of b-values (plus noise), it can be expected that the dynamical behavior is somehow reproduced. More relevant for practical applications is the question how much predictive power the RNN has if it is supplied with only a certain range of parameter values. Is it possible to predict bifurcation points and Lyapunov exponents correctly? To answer this question, we heavily cut the external driving. We cut the probability distribution at b = 0.5: for one RNN we provide only values smaller, for another one larger than 0.5. Then we observe the effect on bifurcation curve and Lyapunov exponent estimate in the respectively "absent" parameter region. The result is shown in Fig. 3 Concluding the results for the Roessler system, we find a comparatively high predictive power of the model obtained. If a whole region of b is not used for training, the result is deteriorated significantly. It is possible to still find a very good prediction of the dynamical regime, if the distance to the values used for training is not too large.

Fitzhugh-Nagumo oscillator -example of an excitable system
In the following we want to study the power of RNN for two systems with different origin and dynamical behavior of oscillations. As a first important class we investigate an excitatory system, namely the Fitzhugh-Nagumo oscillator 43 : where parameters are σ = 0.1, a = 0.7, b = 0.8 and I 0 = 0.25. For the input we use, as above p(t) = q(t), as described by Eq. (3), with = 0.05 and τ = 25.0. The RNN is trained on timeseries comprising of 1000 spikes, over 500 training epochs. The time resolution is ∆ t = 1. For the excitatory oscillations we validate the model on a spike train, given as novel input realization, cf. Fig. 4a. For system (5) the spiking frequency is an important characteristics, in fact relevant for most of the excitatory biological systems modeled. Therefore, we estimate the spiking frequency with respect to the input current I 0 and compare it to the true one, see Fig. 4b. When I 0 is increased, a bifurcation occurs, since above threshold the system "fires" as one says in the context of neural systems. The corresponding first-order phase transition is clearly inferred, with the critical value of the input accurately predicted up to the order 10 −3 , in addition, the estimated frequency values match closely the true ones with mean deviation 10 −2 .

Mackey-Glass equation -example of a delay system
Eventually, we briefly report on the RNN results for the Mackey-Glass system, as a representative for a time-delay system. The Mackey-Glass equation 44 reads: where x θ = x(t − θ), a = 2, b = 1, n = 8 and the time-delay θ = 2. We use input p(t) = q(t), Eq. (3), with = 0.005 and τ = 1.0. In this parameter regime the equation yields a stable limit cycle with a period-2 orbit, see Fig. 5. The length of time series corresponds to 5000 natural periods used as test data over 500 epochs in batches of 100 time points with resolution ∆ t = 0.15. Actually, a delay equation seems to be suited perfectly for a RNN since the delay is inbuilt into the model. Another strong point for the RNN is the rational function which forms the right-hand side of Eq. (6) -as an additive model it is a weakly converging series such that either some network structure or a more general non-additive model is needed. In Fig. 5 we find these expectations confirmed. The dynamics is well reproduced with a mean deviation of 5×10 −2 . We conclude that as well for this important model class RNNs work well, and the strategy to add noise in order to test the phases space is successful, here too.

C. Amount of data and noise study
Any good statistics-based study includes a section on the dependence of the result on the amount of data provided and the sensitivity to noise -we do so in the following paragraphs. We present only results for the Roessler oscillator, Eq. (2). We train in independent runs RNN models with different lengths of input timeseries. Since we already sample the periods fine enough, we vary the number of periods supplied to the RNN in the following way: we keep the product of the timeseries length and number of epochs constant, thereby introducing the same number of data points to the network (500000), i.e. we change the number of occurrences of equal points. We test the range from 15 to 1000 periods and measure the error of the PRC and signal, see Fig. 6. For each set of parameters 100 models are trained and evaluated. The PRC error is evaluated as the L 2 norm of the difference between the true and the reconstructed curve: and similarly for the signal error: where we further have to determine over what interval we evaluate it, ∆.
The results of this evaluation are shown in Fig. 6. Here, we demonstrate the difference between LSTM and GRU cell types, in addition, to underline our previous remark on the poor results for GRU. We use two evaluation criteria for the model quality: the least-squares error of the PRC and the least-squares error of the signal with the respective estimate. After approximately 100 periods, the error in PRC has converged to ca. 0.01. For the least squares deviation of estimate and signal, however, one cannot detect a convergence, there might be a lower error with increasing number The bifurcation regime is very well recovered in that region, beyond b = 0.5 the behavior is qualitatively met, however the exact curve is slightly different. Right: only values b > 0.5 are used for training. The reconstruction for the region b < 0.5 has now significant missing values. I.e. if we do not provide signals that are chaoti enough, a quantitative reconstruction of the bifurcation diagram or the Lyapunov exponent is not possible. However, we want to underline, that for the Lyapunov exponent we still have a surprisingly well qualitative coincidence except for the period-3 window which is not found at all. of periods. It is hard to say at this stage of the research if this is relevant or may be overfitting already. In any case, the GRU cells perform poorly and we can discard them for our questions under investigation. Now for the robustness of the inference against measurement noise. We only present a basic study where we consider the Roessler system, Eq. (2) with input p(t) = q(t), Eq. (3). We fix b = 0.6 and consequently the system to a chaotic state. Then, we add to each timepoint a random uncorrelated Gaussian number with mean 0 and standard deviation 1 to represent strong measurement noise, and train the network on them. We introduce 10000 average periods worth of training data over 500 epochs. The network effectively extracts the relevant dynamics and reproduces the attractor well, see Fig. 7.
We don not systematically study a dependence on noise amplitude, that may be subject to further investigations on quantitative studies on the dependence on provided parameter ranges of b and more details on the other types of system.

III. METHODS
In this section we specify the methods we used to evaluate the properties of oscillatory systems. For each property we write how we computed it from the equations as well as how we computed it from the RNN.

A. Natural period estimation
The period is measured as the time between two successive signal-threshold crossings from bellow when the system is unperturbed. From equations, the time of crossing is accurately estimated using the Hénon trick 45 . When estimating from a network, a linear interpolation from a point before and after the threshold crossing is used.

B. PRC estimation
Firstly the natural period T 0 has to be accurately estimated, see section III A. Then the system in question is weakly and instantaneously perturbed at particular phases Both have the same input realization, depicted with gray (scaled and shifted for being visually comparable to the signal). b: The true spiking rate of system (5) in red and the one inferred from the RNN in green.
FIG. 5. a: Delay embedded trajectory of system (6) in red, its two dimensional projections in purple. For the chosen set of parameters the system has a stable limit cycle, the variation is due to the input. b: The signal generated with Eqs. (2) in red and the RNN reproduced signal in green. Both have the same input realization, depicted with gray (scaled for being visually comparable to the signal). ϕ * , i.e. at times t * = ϕ * 2π T 0 after the begining of a period, p(t) = δ(t − t * ). Then the evoked phase shift is evaluated as where T 1 is the period in which the perturbation arrives and T 2 , T 3 , T 4 , ... the periods that follow. n counts how many periods we wait to evaluate the shift and since we are looking for the asymptotic shift n should be big enough that the PRC does not depend on it, in this paper we used n = 5.
In the case of the network, the time for inputting perturbations is discrete and the best we can do is input perturbation /∆t where ∆t is the time increment between two consecutive points in the unrolled RNN.

C. Maximal Lyapunov exponent estimation
For computing the exponents from the true system we use the standard technique, since we have the dynamical equations.
To estimate the exponent from the RNN a different approach is needed. Suppose we have access to all the variables of the original system x = (x 1 , x 1 , ..., x nx ). In such case the intuitive method can be used: 1. simulate a trajectory x for a long time so it settles to the attractor, 2. start a new trajectory x † = x + p with a small arbitrary perturbation p = δx and evolve both for a short time δt, 3. evaluate the deviation ∆ = x † − x , 4. renormalize the second trajectory for the deviation to have the same amplitude as the one we started with x † = x + δx * ( x † − x)/ x † − x , but keep the direction of the perturbation the same so that the maximal exponents takes over in the course of several repetitions, 5. loop to step 3 and average the quantity 1 δt log(∆/δx) which tends towards the maximal Lyapunov exponent.
Here · stands for the L 2 norm: v = The more general approach concerns cases where we do not have access to all the variables but only a few, in the extreme case only one x 1 -common when dealing with real data. In such case the state of the system has to be characterized with several historical values, w = (x 1 (t), x 1 (t − ∆t), x 1 (t − 2∆t), ...), and then the algorithm above can be used as before. This is the case in section II B 1.

IV. DISCUSSION
The aim of this study was to test the predictive capacity of recurrent neural network applied to different oscillatory systems. One problem common to all oscillators is that the state space collapses to a low-dimensional manifold and therefore any reconstruction only allows the prediction on that inertial manifold. However, if perturbed we can achieve a much better understanding of the system around its attractor. We even can follow and predict a bifurcation outside the range of values which were provided by the data. This is a notable fact and it may as well work for other methods, like symbolic regression.
We have applied the method to a range of oscillatory systems, from a time-delay oscillator with a period-2 orbit (Sec. II B 3, Fig. 5), to an excitatory system (Sec. II B 2, Fig. 4), and finally a chaotic attractor (Sec. II B 1, Fig. 2). We demonstrate that the trained neural networks can be probed for dynamical responses. As typical characteristics of oscillatory systems we estimated the phase response curve 26 (PRC), the spiking rate, and the maximal Lyapunov exponent 27 . Other quantities, such as the Floquet exponent 46 , the amplitude response, the isochronal structure 47 , synchronization properties 48 etc. could be estimated in a similar way. We can say that RNNs provide an effective way of estimating oscillatory properties from timeseries, cf. [40][41][42]49,50 . Our way of applying them to data is novel and should be explored further not only in the context of oscillations. It is, for example, not clear how well RNNs perform for scaling systems like turbulence.
Since each machine learning method depends on data, we performed a statistical analysis on how the size of training data set influences the inference. The training data required for an effective inference proved to be reasonably small, with only a few 10 periods sufficing for reliably estimating the mentioned dynamical systems quantities. We used two popular recurrent network cells in our study: the long shortterm memory cell 33 and the gated recurrent unit 34 . The latter proved to be inferior in performing these tasks. We also tested the inference with the addition of measurement noise and it proved to be robust, see section II C.
Along with this publication, we (RC) published a Python software package, OscillatorSnap 51 , available on the Python Package Index (PyPI) as: oscillator snap. It contains most of the examples shown here as well as an array of high level functions for analysing oscillatory systems, such as, a function that computes the phase response curve or the maximal Lyapunov exponent from dynamical equations as well as from a trained RNN model.