19 June 2000 Physics Letters A 271 Ž2000. 217–222 www.elsevier.nlrlocaterpla Extracting model equations from experimental data R. Friedrich a,) , S. Siegert a , J. Peinke b, St. Luck ¨ b, M. Siefert b, M. Lindemann c,d , J. Raethjen d , G. Deuschl d , G. Pfister c a c Institute for Theoretical Physics, UniÕersity of Stuttgart, D-70550 Stuttgart, Germany b Department of Physics, UniÕersity of Oldenburg, D-26111 Oldenburg, Germany Institute for Experimental and Applied Physics, UniÕersity of Kiel, D-24098 Kiel, Germany d Department of Neurology, UniÕersity of Kiel, D-24105 Kiel, Germany Received 18 April 2000; accepted 2 May 2000 Communicated by J. Flouquet Abstract This letter wants to present a general data-driven method for formulating suitable model equations for nonlinear complex systems. The method is validated in a quantitative way by its application to experimentally found data of a chaotic electric circuit. Furthermore, the results of an analysis of tremor data from patients suffering from Parkinson’s disease, from essential tremor, and from normal subjects with physiological tremor are presented, discussed and compared. They allow a distinction between the different forms of tremor. q 2000 Elsevier Science B.V. All rights reserved. PACS: 05.10.Gg; 05.45.-a; 05.40.Ca 1. Introduction A basic aim of scientific research is to set up reasonable models for considered systems. A suitable model should reproduce the observed quantities and help to gain a deeper understanding of the system. Usually, collected data and known properties of the system, as symmetry relations for example, serve as a basis for the modelling. In contrast to this a general data-driÕen way for formulating suitable model equations for nonlinear complex systems is presented in the following. ) Corresponding author. Fax: q49-711-685-5098. E-mail address: fiddi@theo3.physik.uni-stuttgart.de ŽR. Friedrich.. An important and wide class of dynamic systems can be described by the following differential equation d dt ^ ` _ X Ž t . s g Ž X Ž t . ,t . deterministic part ^ ` _ q h Ž X Ž t . ,t . G Ž t . , stochastic part Ž 1. the Langevin equation w1,2x. X Ž t . denotes the time dependent d-dimensional stochastic vector which characterises the system completely. The time derivative of X Ž t . can be expressed as sum of a deterministic part g and a stochastic part h P G , 0375-9601r00r$ - see front matter q 2000 Elsevier Science B.V. All rights reserved. PII: S 0 3 7 5 - 9 6 0 1 Ž 0 0 . 0 0 3 3 4 - 0 218 R. Friedrich et al.r Physics Letters A 271 (2000) 217–222 where G Ž t . stands for terms of d-correlated Gaussian white noise and where the d = d-matrix h fixes the dynamic influence of the noise on the system. For the functionals g and h no further assumptions have to be made; g can be nonlinear, and therefore also deterministic chaos can be formulated by a Langevin Eq. Ž1.. The investigation of complex systems has shown the necessity to find a description of a system by Langevin equations directly from measured data sets w3x. For this wide class of dynamic systems Ž1. a general method for finding the deterministic and stochastic laws solely by data analysis will be presented in this letter. If this only condition, the describability of the system’s dynamics by an evolution equation like Ž1., is given, no further assumptions or pre-knowledge have to be included in the following analysis. Deterministic and noisy parts of the dynamics can be separated and quantified, and model equations for the dynamics can be set up by the datadriven method. 2. Numerical method The considered class of dynamic systems Ž1. is characterised by the missing of memory effects, i.e. for the time development of the system we need to know only the state of one vector X Ž t . at a given time t and not its evolution in the past. A system with this quality is called a Markovian system w4x. The conditional probability density pŽ x nq 1 ,t q t N x n ,t; x ny1 ,t y t ; . . . . describes the probability of states x nq 1 of the system’s variable X at time t q t under the condition that the system is in state x n at time t, has been in state x ny 1 at time t y t and so on. The Markovian property of a system can be tested by comparing the multiple conditional probability densities with the one-step conditional probability densities pŽ x nq 1 ,t q t N x n ,t .. If both expressions agree the time development of the probability density depends only on the present state and not on its evolution in the past. The assumed qualities of the driving noise terms G as being Gaussian white noise functions can be validated, as well, by looking at the conditional probability density distributions. The central ideas of the presented algorithm will be explained in the following: First, stationary dynamics shall be assumed, i.e. the deterministic and stochastic parts g and h are not explicitly time dependent 1. Every time t i , the system’s trajectory meets an arbitrary but fixed point x in state space, the localisation of the trajectory at time t i q t is determined by the deterministic function g Ž x ., which is constant for fixed x, and by the stochastic function hŽ x . G Ž t i . with constant h for fixed x and Gaussian distributed white noise G Ž t .. With this background in mind, the following relationships which have been proved in a strict mathematical way w5,6x using Ito’s ˆ definitions for stochastic integrals w7x, become understandable: t 1 ™0 t g Ž x . s lim ² X Ž t q t . y x :< X Ž t .sx t 1 ™0 t h Ž x . P hT Ž x . s lim Ž 2. ²Ž X Ž t q t . y x . T = Ž X Ž t q t . y x . :< X Ž t .sx Ž 3. Under the condition that the system’s trajectory meets the point x at time t, i.e. X Ž t . s x, the deterministic part g of the dynamics can be evaluated for small t by the difference of the system’s state at time t q t and the state at time t, averaged over an ensemble, or in the regarded stationary case, averaged over all t s t i of the whole time series with X Ž t i . s x. The limit t 0 can be reached by extrapolation. In a similar way, the stochastic influences can be determined, regarding quadratic terms in the averages w8x. For every point x in state space, that is visited statistically often by the trajectory, deterministic and stochastic parts of the dynamics can be estimated numerically. A schematic illustration of this procedure is shown in Fig. 1. As final step, analytic functions can be fitted to the numerically determined values of g Ž x . and hŽ x . ™ 1 Non-stationary dynamic systems, i.e. with functionals g and h being explicitly time dependent, can be investigated by a moving window technique. An analysis window is drawn along the whole time series in overlapping steps. In each window, the system’s dynamics is assumed to be stationary. When interpreting the results, the quality of the algorithm of calculating weighted averages of the real functionals has to be taken in mind. R. Friedrich et al.r Physics Letters A 271 (2000) 217–222 219 scribed above. Afterwards, the numerically determined results and the expected results according to the system’s evolution equations are compared. The electric circuit is shown as schematic illustration in Fig. 2. Its dynamic equations are given by the following evolution equations, where the deterministic part is known as Shinriki oscillator w9x: 1 1 f Ž X1 y X 2 . Ẋ 1 s y y X1 y R NIC C1 R 1 C1 C1 1 q G Ž t. R NIC C1 s g 1 Ž X 1 , X 2 . q h1 G Ž t . Ž 4. f Ž X1 y X 2 . 1 Ẋ 2 s y X C2 R 3 C2 3 ž Fig. 1. Schematic illustration of the presented algorithm for analysing noisy data sets and calculating the deterministic and stochastic parts of the underlying dynamics. In the upper part of the figure a trajectory of a chaotic oscillator ŽShinriki. with dynamical noise is shown. From the trajectory all parts are selected which run exactly or within some limits at any time t i through a certain point x in state space, i.e. X Ž t i . s x q D x. The distribution of the values x˜ taken at the next step X Ž t i qt . of these trajectory parts is given by a Gaussian function with mean x q g Ž x .t and standard deviation hŽ x .'t , where g and h are defined by Eq. Ž1.. g Ž x . is the mean velocity of the trajectory at point x, hŽ x . represents the fluctuations of the velocity at this point. The conditional probability density pŽ x,t ˜ qt N x,t . can be evaluated by this way. / s g 2 Ž X1 , X 2 , X 3 . Ž 5. R3 Ẋ 3 s y Ž X2 y X3 . L s g 3 Ž X2 , X3 . . Ž 6. X 1 , X 2 and X 3 denote voltage terms, R i are values of resistors, L i and Ci stand for inductivity and capacity values. The function f Ž X 1 y X 2 . denotes the characteristic of the nonlinear element. The quantities X i , characterising the stochastic variable of the Shinriki oscillator with dynamical noise, were measured by means of a 12 bit ArD converter. Our analysis is based on the measurement of 100.000 data points, the affiliated trajectory is shown in the upper part of Fig. 1. The measured three-dimensional time series were analysed according to the algorithm described above. The determined deterministic dynamics – expressed in order to formulate model equations for the investigated system. 3. Analysis of electronic data Next, the application of the method to data sets from two experimental systems will be presented. As first example, a chaotic electric circuit has been chosen. Its dynamics is formed by a damped oscillator with nonlinear energy support and additional dynamic noise terms. In this case, well defined electric quantities are measured for which the dynamic equations are known. The measured time series are analysed according to the numerical algorithm de- Fig. 2. Illustration of the investigated electric circuit. The structure of the dynamics is a chaotic oscillator with dynamical noise. 220 R. Friedrich et al.r Physics Letters A 271 (2000) 217–222 g 1Ž x 1 , x 2 s 0. is drawn in part Žb.. Additionally to the numerically determined results, found by data analysis, the expected vector field and curve ŽEqs. Ž4., Ž5.. are shown for comparison. A good agreement can be recognized. 4. Analysis of tremor time series Fig. 3. Cuts of the function g Ž x . reconstructed from experimental data of the electric circuit Žillustrated in Fig. 2. in comparison with the expected functions according to the known differential Eqs. Ž4. – Ž6.. In part Ža. the cut Ž g 1Ž x 1 , x 2 ., g 2 Ž x 1 , x 2 , x 3 s 0.. is shown as two-dimensional vector field. Thick arrows represent values determined by data analysis, thin arrows represent the theoretically expected values. In areas of the state space where the trajectory did not go during the measurement no estimated values for the functions exist. Figure Žb. shows the one dimensional cut g 1Ž x 1 , x 2 s 0.. Points represent values estimated numerically by data analysis. Additionally, the affiliated theoretically curve is printed, as well. by the deterministic part of the evolution equations – corresponds to a vector field in the three-dimensional state space. For presentation of the results, cuts of lower dimension have been generated. Part Ža. of Fig. 3 illustrates the vector field Ž g 1Ž x 1 , x 2 . , g 2 Ž x 1 , x 2 , x 3 s 0. . of the reconstructed deterministic parts affiliated with Ž4., Ž5.. Furthermore, the one-dimensional curve In the second example, data originating from neurophysiology are investigated. Time series of hand tremor have been recorded in normal subjects with physiological tremor ŽPT., in patients with essential tremor ŽET. and in patients suffering from Parkinson’s disease ŽPD.. The data sets were measured by a lightweight piezoresistive accelerometer with a sampling rate of 800 Hz. The outstretched hand was supported at the wrist. For the analysis 24 000 data points per experiment were used. The data sets were scaled to variance 1. It is assumed that the underlying dynamics of the measured time series X 1Ž t . can be expressed by a two-dimensional Langevin equation Ž1.. Using a time delay method a second time series X 2 Ž t . was created out of the measured time series X 1Ž t .. Afterwards, the algorithm described above was applied. Up to now 20 data sets have been investigated, measured at patients whose disease can be clearly classified within one of the three groups: physiological tremor, essential tremor and Parkinson’s disease. Further data sets will be investigated, but already now a clear structure of the dynamic behaviour appears. Typical results of the analysis for the three tremor groups together with extracts of the investigated time series can be seen in Fig. 4. Whereas the deterministic parts of the dynamics in the case of physiological tremor can be described as a fixed point with large damping and weak rotation, the deterministic parts of the dynamics in the case of essential tremor are characterised by a fixed Fig. 4. Typical results of an analysis of tremor data achieved by the presented algorithm. In part Ža. variance normalised accelerometer data from a normal subject with physiological tremor was investigated, in part Žb. from a patient with essential tremor and in part Žc. from a patient suffering from Parkinson’s disease. In every part a subsequence of the used time series and the numerical results for the deterministic parts of the dynamics, drawn as vector diagram, are shown. The curves are trajectories integrated along these vector fields. Whereas in part Ža. and Žb. the underlying deterministic parts of the total dynamics can be described as fixed points the deterministic parts of the dynamics in the case of Parkinson’s disease are characterised as periodic limit cycle. The dynamics of physiological and essential tremor can be separated by their damping, i.e. the number of cyclings of the trajectory on its way to the fixed point. R. Friedrich et al.r Physics Letters A 271 (2000) 217–222 221 222 R. Friedrich et al.r Physics Letters A 271 (2000) 217–222 point, as well, but with much weaker damping and larger rotation resulting in a lower inwards spiralling motion. In the case of Parkinson’s disease in contrast, the deterministic parts of the dynamics are governed by limit cycles. The result of the analysis concerning the physiological tremor is in accordance with earlier studies on model equations for the dynamics of physiological tremor. In Ref. w10x it was shown that these dynamics can be explained by an ARŽ2. process driven by the uncorrelated firing of motoneurons in the arm muscles. 5. Conclusions In this letter a method is presented which enables one to extract the underlying dynamic equations from observed data. It is possible to estimate the deterministic and the stochastic parts of this dynamics. For validation, the method was applied to data from a chaotic electric circuit, first. The dynamic equations of this physical system were known and allowed a comparison of numerically determined and expected results. Furthermore, the algorithm was applied to tremor data measured in normal subjects with physiological tremor, patients with essential tremor and patients with Parkinson’s disease. The results allow a distinction of the different diseases. As this method allows a better understanding of the underlying dynamics of tremors more insights into their pathophysiology are expected. Concluding, it is expected that the presented method will have a huge field of possible applications in science and might become a new standard tool in the analysis of time series of nonlinear, complex systems. Acknowledgements The authors acknowledge helpful discussions with A. Kittel and Ch. Renner, Department of Physics, University of Oldenburg. References w1x H. Haken, Synergetik, Springer, Berlin, 1990. w2x H. Haken, Advanced Synergetics, Springer, Berlin, 1983. w3x H. Haken, Information and Self-Organization, Springer, Berlin, 1988. w4x P. Hanggi, H. Thomas, Phys. Rep. 88 Ž1982. 207. ¨ w5x A.N. Kolmogorov, Math. Ann. 104 Ž1931. 415. w6x H. Risken, The Fokker-Planck equation, Springer, Berlin, 1984. w7x K. Ito, ˆ Nagoya Mth. J 1 Ž1950. 35. w8x S. Siegert, R. Friedrich, J. Peinke, Phys. Lett. A 243 Ž1998. 275. w9x M. Shinriki, M. Yamamoto, S. Mori, Proc. IEEE 69 Ž1981. 394. w10x J. Timmer, Int. J. Bif. Chaos 8 Ž1998. 1505.