Phase space reconstruction from a biological time series. A PhotoPlethysmoGraphic signal a case study

In the analysis of biological time series, the state space comprises a framework for the study of systems with presumably deterministic properties. However, a physiological experiment typically captures an observable, or, in other words, a series of scalar measurements that characterize the temporal response of the physiological system under study; the dynamic variables that make up the state of the system at any time are not available. Therefore, only from the acquired observations should state vectors reconstructed to emulate the different states of the underlying system. It is what is known as the reconstruction of the state space, called phase space in real-world signals, for now only satisfactorily resolved using the method of delays. Each state vector consists of m components, extracted from successive observations delayed a time t. The morphology of the geometric structure described by the state vectors, as well as their properties, depends on the chosen parameters t and m. The real dynamics of the system under study is subject to the correct determination of the parameters t and m. Only in this way can be deduced characteristics with true physical meaning, revealing aspects that reliably identify the dynamic complexity of the physiological system. The biological signal presented in this work, as a case study, is the PhotoPlethysmoGraphic (PPG) signal. We find that m is five for all the subjects analyzed and that t depends on the time interval in which it evaluates. The H\'enon map and the Lorenz flow are used to facilitate a more intuitive understanding of applied techniques.


I. INTRODUCTION
Dynamic systems, as could be any physiological system, are mathematically characterized by differential equations. The identification of the simplest physiological model, to describe the physiological temporal evolution, requires specifying the minimum dimension of the dynamic system, that is, the number of dynamic variables or differential equations involved in the evolution of the system. The dynamic variables compose the components of a state vector, in the state space 1 , that describes the dynamics of the system. At each time instant, the state vector is in a different position (an isolated point in the state space); the chronological evolution of these points draws up a trajectory in the state space. When the trajectory extends to infinity, it is known as an orbit 2 .
The time evolution of each dynamic variable involves the most natural way to characterize any dynamic system; this is nothing but the usual representation of the value of each dynamic variable as a function of time 3 . Another way to describe a dynamical system graphically is to replace the timeindependent variable, typically, the variable time t, with another dynamic variable of the system 4,5 . In this case, each point on the graph, in a two-dimensional coordinate system, represents the system state in a given time instant; with a third dynamic variable, a three-dimensional coordinate system is available, the highest possible visual capacity. In a hypoa) Electronic mail: javier.depedro@uah.es b) Electronic mail: anamaria.ugena@upm.es c) Electronic mail: anapilar.gonzalez@upm.es thetical abstract space, the coordinate system must cover all dynamic variables of the system under study. Each point of the theoretical graph represents a system state at a given time instant. This graphic construct designates what is called the state space, where each point denotes a state vector. If the dynamic variables that make up the state vectors evolve according to the established physical laws (primitive concepts), as in the real-world signals, see the biological signals, the theoretical and abstract state space becomes a space restricted to states with physical significance; it is from now on the socalled phase space. The graphic display of different trajectories traced by the state vectors in the phase space is termed phase diagram.
The graphical analysis of a dynamic system is reduced to a state space of at most three-dimensions, implicitly including the time variable 6 , without forgetting that other unknown variables, apart from the three considered, can be significant in the system dynamics. Since most real-world physical systems are nonlinear, with an ever-present coupled noise, from a graphic analysis, with the appropriate dynamic variables, at least the deterministic structure of the underlying system dynamics can be deducible, not apparent in a possible erratic evolution in the time domain. A more complex study, although less visual, with many more variables, requires a mathematical formalism not yet fully consolidated 7 .
In physiological terms, typically, in a clinical trial, it is customary to acquire a unique biological signal identifying each physiological response, so that not all the dynamic variables involved in the dynamics of the system under study are available. Thus, each biological signal represents a response of a physiological system; each response distinguishes the time evolution of a dynamic variable. Accordingly, from the scalar measurements that define the acquired biological signal, the different state vectors must be generated in order to recreate the different states that characterize the physiological system dynamics in question. The phase space reconstruction allows to faithfully reproduce the evolution of the system under study from a simple biological signal, in the absence of the dynamic variables involved.
Usually, observations or measures at regular time intervals, i.e., one every ∆t seconds, is an example of what in academic slang know as a time series. The most common method to define each state vector is to use delayed versions of the observations as state vectors components. In the simplest case, in a two-dimensional coordinate system, as illustrates Fig. 1 for the two representative dynamical systems that we use in this paper to explain the methodology employed, the ordinates axis represents the measurement value x n at time n∆t and the abscissas axis the measurement value x n−1 at time (n − 1)∆t. Hence, each point in the plot identifies the x n = (x n−1 , x n ) state vector. At the most general level, the number m of components of the state vector and the delay τ between components is variable, that is, x n = (x n−(m−1)τ , . . . , x n−τ , x n ). From a sequence of scalar measurements, it can reconstruct every state vector, in m-dimensional phase space, following the lag or delay method 8 . The state space reconstruction quality depends on the proper estimation of the parameters τ and m, considering the time delay embedding theorem [9][10][11] , which gives the ideal conditions under which the chaotic dynamical system can reconstruct from a sequence of observations of the state of the dynamical system.
The state or phase space reconstruction, in contrast to the traditional signal analysis in time or frequency domains, constitutes the cornerstone of the time series analysis in terms of nonlinear dynamics 12,13 , and it is the first step in the analysis process of the real functional nature that reflects the physical system in question. In doing so, special attention must be given to determine the optimum values of τ and m 14 efficiently, so that the state space reconstruction faithfully allows to characterize the dynamics of the original system. We detail the most common techniques, elegant in its simplicity, which enable a first approximation of the state space reconstruction and then apply them to a particular biological signal, the Pho-toplethysmographic signal, extracting interesting conclusions that lead to promising future studies.
In a typical trajectory, over the phase space, dynamic variables time evolution can tend to infinity or confined to a bounded region of the phase space. If the dynamic system is dissipative, once the transient response finishes, its dynamics will tend to a subset of the state space called attractor. This subset is invariant under the dynamical system evolution. In a chaotic system, the attractors describe very complex geometric objects, having a typical fractal structure; they are the so-called strange attractors.
The attractor geometry provides valuable information not only on the dynamic nature (Lyapunov exponents) of the underlying physical system but also about the structural complexity sustaining that dynamics (dimensionality), hence the vital importance of a successful phase space reconstruction. The connections between these factors go beyond the scope of this paper, though we will deal properly with these issues in future communications.
In this paper, we outline the due process to be followed in the phase space reconstruction from one scalar time series, applying the methodology to a relatively well-known and easily accessible biological signal, the Photoplethysmographic signal. In section II, we introduce several graphical methods to intuitively assess the approximate determinism present in the dynamic system evolution, making the state space reconstruction meaningful. After, in the next section, in section III, we explain in general terms the main mechanisms underpinning the state space reconstruction, focusing mainly on how to properly determine the reconstruction parameters, the lag τ and the embedding dimension m. In section IV, we briefly describe the basic characteristics of the biological signal used in this work, the PhotoPlethysmoGraphic (PPG) signal, and show the obtained results in line with the methodology referred to above. Finally, in section V, we analyze and interpret the obtained results laying the ground for forthcoming works.

II. DETERMINISM AND STATE SPACE REPRESENTATION
From a linear perspective, the erratic (or irregular) behavior of the one system response is due to a random external component. The chaos theory finds that a random input is not the only cause to get out an irregular output. A simple deterministic equation, as it is with nonlinear chaotic systems, can generate irregular data without the contribution of any external input.
In a deterministic system, according to the deployed dynamical systems theory 15 , once known its current state, the future states are entirely determined. The system dynamical and geometric properties all are to some extent included in the state space representation. Equations of motion that describe the behavior of a physical system as a function of time are deployed in the state space. Henceforth, the drawn trajectories in the state space represent the time evolution that the state vector is undergoing over time. For this reason, the state space or phase space can be used to get an approach to the rules that govern the dynamical system evolution. In an experimental setting, the system equations are not available. Usually, it is had solely with one system response projection in a one-dimensional space, a scalar time series, which represents the acquired measurements of the system in question. The state space reconstruction method aims to reconstruct the state vectors from the time series, so that the time evolution of these vectors replicates a dynamics equivalent to that of the original system, as shown in Fig. 2. A time series reflects roughly the one dynamical variable time evolution, but as J. Doyne Farmer once said the evolution of a variable depends on other variables of the system or, in other words, its value is within the history of other variables with which it interacts 16 . Ideally, a system of first-order ordinary differential equations taking action on the state space 2,12,17 defines a dynamical system. A set of possibly infinite states, and certain transition rules that specify how the system moves from one state to another can describe several systems. For a finitedimensional state space R m , an state is defined by a vector x ∈ R m . Then the dynamics of the system can be described by an m-dimensional map or by a system of m first-order ordinary differential equations called flow. In the first instance, such as the Hénon map, the time is a discrete variable, while in the latter, as is the Lorenz flow, the time is a continuous variable, For each initial condition x 0 , or x(0), the solutions to equations (1) and (2), a sequence of points x n , o x(t), respectively, describe a dynamic system trajectory. With its time evolution a typical trajectory can tend to infinity or confine to a bounded region in the state space. All the initial conditions that lead to the same asymptotic behavior of the observed trajectory is known as basin of attraction 18 .
The trajectory described in the state space, from a single observable, is a presumable indication of the presence of a deterministic behavior if the states not arrange as a point cloud, despite an erratic appearance in the time domain. In a 2D phase diagram the variable time, on the abscissa axis, is replaced by the observable value in a prior time determined by the parameter τ, as discussed in the section I; in a 3D phase diagram except for the first axis, all other axes pick up previous values of the observable variable separated from each other a time τ∆t.

A. Poincaré section
It covers the geometric figure that describes the evolution of a trajectory in a cross-section of an attractor, transversal to the flow or bundle of trajectories, as highlighted in the inner plane in Fig. 3a, gray coloured. Fig. 3b shows the points contained in the Poincaré section. The Poincaré section aims at facilitating a complementary perspective of the general dynamics of the system (internal structure of the attractor), contributing to identify the type of attractor 19 . The Poincaré section describes a pattern of points that can be labeled chronologically, although the time between successive intersections with the section is irrelevant. It must not be confused with a Poincaré map, also called a return map, which reconstructs the temporal sequence of spatially arranged points on the Poincaré section. In a chaotic attractor, the Poincaré section usually describes a complex geometry, as reflected in Fig. 3b, with a very fine distinctive structure, adopting shapes with texture or with multiple layers (fractals).

B. Other plots
The following items summarize other methods of visualizing a possible hidden determinism.
a. Next-amplitude plot. From successive maximums detected in time series, each one M n is represented (abscissa) versus its immediate or subsequent successor τ times M n+τ (ordinate) ahead. A well-defined curve, like the one in Fig. 4a, could reveal the presence of chaos, although the noise could mask a correct interpretation.
b. Difference plot. The graphic's coordinates are delayed differences between successive observations, whether immediate or separated τ number of times. On the abscissa axis it is represented ∆s n = s n+1 − s n and on the ordinate axis the next difference ∆s n+τ = s n+τ+1 −s n+τ . In the simplest form, the first-difference plot, with a delay τ = 1, on the abscissa it is represented s n+1 − s n and on the ordinate s n+2 − s n+1 . The presence of an infinitely continuous curve, as illustrated in Fig. 4b, evidences a high degree of underlying determinism. Hénon next-amplitude plot Hénon first-difference plot

III. PHASE SPACE RECONSTRUCTION METHOD
It was in 1980 that Packard et al. 20 first tackled the problem of how to connect the phase space or state space vector x(t) of the dynamical variables of one physical system to a possible time series {s n } measured in any experiment.
A time series is a sequence of scalar measures s n of an observable, acquired at regular times ∆t. The different measures depending on the current state of the system, where ψ represents an observable measurement function and η n the measurement noise which characterizes the random nature of the imprecision of the measure 21 .
An m-dimensional reconstruction of the s n state vectors, from scalar measurements, is given by The time interval between adjacent coordinates of the state vector is τ∆t, and it is known as lag or delay time. Formally It has been shown that if m is higher than two times the capacity dimension of the attractor (dim C ), a one-to-one correspondence between the reconstructed attractor and the real attractor is guaranteed, no matter how large the dimension of the original state space is 10,11 .

A. Lag or delay time selection
In theory, with an unlimited number of measurements without noise, any delay is equally valid. In practice, with real experimental data, the geometrical and dynamic properties that characterize the attractors differ from each other, depending on the delay value chosen 17 . If the sampling time ∆t is tiny, the values of consecutive samples are very similar, s n ≈ s n+τ , and with τ = 1 the represented points contain much redundant information. Also, with a small delay, the noise level can blur or hide any local geometric structure. Regardless of the form that the attractor can take, the graphic arrangement of perfectly related pairs of points describes a straight line, as an identity function (a diagonal line of 45 • ) without any meaning, as illustrated in Fig. 5a. In order for the distribution of points to provide some geometric significance, the values of the components that make up the coordinate system must be sufficiently independent. In the real state space, in which the coordinates correspond to different variables, this condition is satisfied.
In the same vein, a delay too large can be counterproductive. With a significant delay the dynamic relationship between the values of the variable disappears, so that, as the delay value increases, the geometric structure of the points becomes more complex and diffuse until finally the points are dispersed randomly on the state space (see Fig. 5b). In Fig. 6 the Lorenz flow, with the correct delay time τ = 16, is represented. Next subsections describe how to calculate it.

Autocorrelation coefficients
The autocorrelation coefficients R τ measure the correlation degree of a variable with itself at different instances. Its constituents are autocovariance and variance.
where N represents the time series length, and s the arithmetic mean of all the time series observations. According to the equation (5), with a minimum delay τ = 0, the coordinates of each point are identical (s n = s n+τ ), and the autocorrelation is maximum R τ=0 = 1. As delay increases, autocorrelation decreases until it eventually is reduced to zero, or, as happens in real data, there is a fluctuation around zero, within a narrow margin, due to the noise in the data. In short, the autocorrelation coefficients range from +1 y −1. Maximum values, R τ = ±1, are perfectly correlated data; minimum values, R τ = 0, correspond to uncorrelated data.
Autocorrelation coefficients, for successive delays, conform the autocorrelation function. It is also known as the spectral autocorrelation coefficient 22 (it is recommended to use N/4 delays with a time series of N > 50 observations). The graphical representation of R τ is known as correlogram. To a certain extent, a correlogram shows the type of regularity in the data. In the case of trendless and uncorrelated data, 95% of the autocorrelation coefficients are contained, in theory, in a band of ±2/ √ N, around zero. About 5% of them can exceed the indicated limit without lost its condition of uncorrelated data 23 . From the autocorrelation function, it is possible to adopt some criteria for selecting the optimal delay. Among the various selection criteria, the most commonly used determines as optimal delay the first zero crossing of the autocorrelation function, as illustrated in Fig. 7a. As an example of an alternative approach, in climate-related attractors 24 , other criteria apply the one known as correlation time (when R τ decreases to a value of 1/e ≈ 0.37).

Mutual information
Autocorrelation coefficients consider the degree of the mutual relationship on a linear basis, inappropriate in nonlinear systems 14 . In these situations, it is best to use the mutual information, as it enables to determine, on a probabilistic basis, to what extent two values of the same variable, measured at different time instants, relate to each other. For example, if the point coordinates are identical, with τ = 0, these represent the same information and, therefore, one coordinate accurately predicts the other, or, otherwise, the amount of information that one coordinate contains about the other, mutual information, is maximal.
Whenever there is any relationship between two values of the same variable, acquired at different time instants, one value contains information about the other value. That is, one value helps predict the other value or, in other words, the knowledge of one value reduces the uncertainty of the other value. The uncertainty reduction is called mutual information. If X denotes a time series, {s n }, ergo, s 1 , s 2 , . . . , s N , with N observations, and Y the delayed time series, {s n+τ }, ergo, s 1+τ , s 2+τ , . . . , s N , con N − τ observations, where τ indicates the delay, the average mutual information I Y;X between both time series can be expressed, in probabilistic terms, as where N c is the number of cells containing points, with nonzero probability, and N r is the number of routes, s i s i+τ , in the state space. Equation (6), in terms of entropy, rewrites as where H X is the entropy of X, H Y the entropy of Y, and H X,Y the joint entropy of X and Y. So somehow, the mutual information involves a measure of predictability of the system, that is, a measure of the degree of knowledge of s n+τ noted s n . In order to reconstruct an attractor, as faithfully as possible, minimum mutual information and delay value are required. As with the autocorrelation function, the mutual information decreases, as the delay increases, until it eventually becomes zero. The first minimum of the mutual information includes a possible criterion for selecting the optima delay value 8 , as shown in Fig. 7b. One limitation of this method is that many points are needed if we are to get a consistent result.

B. Embedding dimension selection
There is no rule of thumb to set the minimum reconstruction dimension m, and none of the published proposals is widely accepted, so it is a good idea to use more than one method on the same data. Among all possible techniques, the correlation dimension and false nearest neighbors stand out, though principal components analysis, as a preliminary attempt, may shed some light.
From a scalar time series, {s n } N n=1 acquired, and with the proposed τ value as stated above, N − (m − 1)τ vectors with m component per vector are defined, where each component symbolizes an alleged dynamic variable of the physical system. To a certain extent, we initiate a multivariate analysis with m time series data obtained by (4).

Principal Component Analysis (PCA)
A first approach to the estimation of the parameter m comes from the Principal Component Analysis (PCA), an application of linear algebra. Fig. 8a shows the results for Lorenz attractor. This method can help to extract relevant information from seemingly complex data, a priori reducing the number of dimensions needed to characterize the dynamics present in the data, that is, to identify the basis that best represents a noisy data set, detecting those redundant dimensions that record the same dynamic information 13 . This basis, to a certain extent, filters the noise and improves the recovery of the hidden dynamics in data. In all, this linear composition aims to find the smallest subspace (hyperplane) that contains roughly the attractor.
The PCA method takes up an initial reconstruction dimension m much higher than the one that supposedly characterizes the attractor. The experimenter intuition plays an important role in deciding on a greater or lesser initial value of m. It points the directions of the reconstruction space that show more significant variations (variances) in the data, discarding those other directions in which small variations may come from fluctuations of the coupled noise on data. Discarding those less important directions, apart from reducing the effect of noise, particularly white noise, makes it possible for appearance of a more simplified dynamics from a very high dimensional space. Hence, it enables a dimensional reduction. The main disadvantage of this method is subjective nature when determining the value of m. With real data, or even with an unfortunate initial choice of τ and m, the differences between the variances of the different dimensions are not always so evident, and it is necessary to define a somewhat arbitrary threshold that classifies the dimensions as primary or secondary. Furthermore, there is no guarantee that the reconstruction will always be optimal, since the reconstruction method relies on a non-parametric analysis, and sometimes it is unable to distinguish between a chaotic signal and the noise itself when they have a similar power spectrum 25 .

Correlation dimension
The correlation dimension involves the most usual measure of dimension, mainly because of its computational efficiency. It is based on spatial correlation 26 . The general procedure includes point counting within a distance ε, evaluated for each of the points that make up the state space. The normalized total number of points, for a particular distance ε, is called the correlation sum, the estimator of the correlation integral, where N is the time series length. For large enough N, which applies for a N value of several hundred data; according to the mathematical formalism, where Θ is the Heaviside step function, so, Θ(x) = 1, for x ≥ 0, and Θ(x) = 0, for x < 0, and · symbolizes the maximum norm of one vector. The equation (10) is computed for different values of ε, with the vectors obtained from the equation (4) for each particular case of m. If the data obtained, C ε versus ε, are arranged in a straight line in a log-log plot, especially for the intermediate values of ε, that means that the correlation sum follows a power law. For extreme values of ε, the statistical estimation of C ε is not reliable, and, therefore, the points obtained usually move away from the straight line. Consequently, a straight line, in a log-log plot and for a central region of values of ε, suggests a power law, which, in this case, obeys where the exponent dim C ε ≡ D 2 (ε), often a non-integer number, denotes the slope of the line, and it is called correlation dimension.
For different values of m are generally obtained different values of the slope of the line. The minimum correlation dimension value, for which additional increments of m do not clearly modify its value, implicitly defines the appropriate reconstruction dimension, say, m > 2 D 2 , in accordance with Takens's theorem 10 , so that the attractor can fully deploy its dynamics (see Fig. 8b). If the correlation dimension grows continuously with each value of m initially chosen, the data suggest a random behavior.

False nearest neighbors
In a reconstruction space of very low dimension, two points appear to be closer to each other than they are. Two points are considered as real neighbors if the distance between them remains constant as the reconstruction dimension increases. Conversely, the distance between false neighbors continues to increase as long as the reconstruction dimension remains too low.
The underlying principle behind the method is to look for points in the time series that are neighbors in the reconstruction space, but that should not be, as their future time evolution is very different 27,28 . The distances between all points, for different consecutive reconstruction space dimensions, mdimensional and (m + 1)-dimensional spaces, must be estimated. If the ratio between both distances is higher than a threshold r, it says that the neighbors are false neighbors.
If the standard deviation of the data is σ , and it uses the maximum norm, for computation speed reasons, the percentage of false neighbors χ amounts to where s (m) k(n) is the nearest neighbor to s n in the m-dimensional reconstruction space. The subscript k(n) indexes the time series element, with k(n) = n, for which s n − s k(n) is minimal.
The first term of the numerator in equation (12), within the summation, is equal to one, if the nearest neighbor is false, that is, if the distance increases by a factor greater than r when the reconstruction dimension is increased by one, from m to m + 1. The second term drops those pairs of points whose initial distance is already greater than σ /r, since, by definition, they can not be false neighbors, as, on average, the points cannot be further away than σ . Therefore, in the calculation, these points should not be taken and, thus, appear in the denominator's normalization factor. The right reconstruction dimension is that for which the percentage of false neighbors falls to approximately zero, as shown in Fig. 8c. Once it reaches the dimension, it assumed that the attractor embraces its true spatial configuration, at least from a topological perspective. The results may depend on the chosen delay τ. In any case, if it seeks to be able to differentiate the chaos from the noise, it is essential to validate the results using surrogate data testing 29 .

IV. BIOLOGICAL DATA
The health monitoring by non-invasive means has attracted the attention of medical specialists for some years now since it allows to advance preliminary diagnoses on possible pathological dysfunctions, whether chronic or transient. Various scientific evidences sustain the relationship between the biological signals generated by the human body and the health status of the individual, as skin temperature (ST), electrodermal activity (EDA), pulse wave or photoplethysmography (PPG), electrocardiography (ECG), electromyography (EMG), respiration (Resp.), pupil diameter (PD), electroencephalography (EEG), blood pressure, among others. The information provided by all these biological signals depends on the knowledge available about the underlying physiological processes.
The study of the dynamics of all these biological signals will help to understand better the physiological system that generates them. Also, it would help to get an insight about how different dynamic variables couple, as, for instance, the heart and respiratory frequencies, in order to keep the physiological system in a perfect condition, preventing its deterioration towards a pathological state. In this paper, we focus only on a single biological signal, the photoplethysmographic (PPG) signal; forward publications will describe the results for more biological signals.

A. Results on the PPG signal case study
We have chosen the PPG signal because it is easily accessible. A pulse oximeter consists of a light emitter and a photodetector that collects and records (pulse or PPG signal) the loss-scattering and absorption-that a beam of light undergoes when it passes through or is reflected, a human tissue 30 . It allows detecting blood volume changes in the microvascular bed of tissue-in our case, the middle finger of the left hand-, obtaining valuable information about the cardiovascular system and, on the whole, about the cardiorespiratory system. Given the simplicity of its non-invasive ac-commodation, in addition to its low cost, a pulse oximeter is very useful in biomedical applications for clinical 31 and sports environments 32 . As with other biological signals, characteristics extracted of the PPG signal allow that to no small extent identify ideal health conditions and their possible deviations. Thus, indicators associated with different pathologies could be established, which anticipate its severity according to the causes that gave rise them. Typically, these indicators in the case of PPG signal were based on the morphology of the signal rather than on its dynamics 31,33 ; we think that by studying dynamic aspects of the PPG signal, the physiological system that generates it can be better understood.
With the aim of examining the PPG signal dynamics, we use the PPG signals (to show in this paper only five individuals chosen randomly) from a total of 40 students, between 18 and 30 years old and a non-regular consumers of psychotropic substances, alcohol or tobacco, selected to participate in a national research study 34,35 ; the five PPG signals show in this paper are the same as those used in previous research that confirms its predominantly quasi-periodic behavior for small timescales in healthy young people with the 0-1 test 36 . To show more than five subjects will not clarify the proposed method. The results are similar. Remember that the fundamental frequency of the PPG signal is typically around 1 Hz, depending on heart rate (0.5-4 Hz, first and second harmonic), and respiratory activity at 0.2-0.35 Hz. We apply a Butterworth bandpass filter with cutoff frequencies at 0.1 and 8 Hz, in order to avoid high-frequency noise and to some extent, motion artifacts 37 . All signals were captured from the middle finger of the left hand and sampled at a frequency of 250 Hz 34 , say, sampling time ∆t = 4 ms.
The study methodology of the PPG signal to finally get the phase space reconstruction and draw interesting conclusions about the underlying physiological dynamics follow the next steps: first, we take 15,000 points which correspond to one minute recorded signal and study its graphical representation, considering different options for different graphic representations, as shown in Fig. 9 (in the accompanying video clips we can watch the evolution of the attractor geometry for different lags, for both 2D and 3D phase diagrams), to sustain the presence of a latent dynamic structure. Second, we try to figure out the attractor more akin to the original, based on the optimal lag τ, see subsection III A, according to the Autocorrelation Function (AF) (Fig. 10a) and the Mutual Information (MI) (Fig. 10b), and on the optimal embedding dimension m, see subsection III B, in accordance with the Principal Component Analysis (PCA) (Fig. 11a), the Correlation Dimension (D2) (Fig. 11b), and the False Nearest Neighbors method (FNN) (Fig. 11c). We apply the methodology to study 10 minutes of the PPG signal from each subject, and the results are similar to the ones shown in this paper.
We have seen in the Lorenz flow, according to the AF, that the criterion of 1/e is approximately similar to the MI result because the attractor is chaotic and both linear and nonlinear correlations decay rapidly, as illustrated in Fig. 7. In the PPG signal, the criterion of the first zero crossing is closer to the MI result, see Fig. 10, because the PPG signal is predominantly quasi-periodic on small timescales 36  stronger. Either way, the criterion of the first minimum of MI is adopted because it ensures statistical independence between the values, in linear and nonlinear terms. The parameters τ and m calculate for five subjects, with results from each of the methods previously described, are in Tab. I. The optimal τ depends on the subject considered and on the time interval in which it estimates. Probably the range of optimal τ values varies with each subject and could be a hallmark of each individual. It should make a careful study with more subjects.
Once we have the parameter τ for each subject, calculated by the MI criterion, we now proceed to the calculation of parameter m. Fig. 11 visually represent the results for a sample PPG signal, which corresponds to subject number 4 of the Tab. I. We see that multivariate analysis, such as the PCA linear method, does not provide the ideal value of m, although it is not unreasonable, mostly because physical systems do not usually have very high dimensions. Correlation dimension gives the value of m > 2 D 2 , that, except for the first subject, amounts to a value of m = 5. For the FNN method, all calculations agree with m = 5. This later method is a good estimator of m because in resorting to topological aspects of the reconstructed state space, the effect of noise on algorithmic computation is minimal. We conclude that all individuals have the same embedding dimension m = 5, although we think that it may vary depending on the subject's physical and psychological states. Consequently, a system of five first-order differential equations describes the dynamic system that generates a typical PPG signal from a healthy young individual, in which various parameters could contribute to the mechanism of autoregulation of the underlying physiological process. Fig. 12 shows how the attractors of the subjects are different for the calculated values of τ. 38 Ongoing studies on the effect of changes in the parameters τ and m on the dynamics of the underlying physiological process will shed any light in this regard.

V. CONCLUSION
The well-known phase space reconstruction serves as a starting point for the analysis or modeling of the dynamics of physiological systems reconstructed from a single biological signal. The dynamic that describes a photoplethysmographic signal is deterministic. The embedding dimension m, for healthy young subjects, is five, regardless of the subject and the time interval estimation. The optimal lag or delay time τ opt depends on the subject and calculated time interval. A further study, with individuals of different ages and with different proven pathologies, will be carried out to confirm if the embedding dimension remains the same. In addition, it will be thoroughly examined if the variation of τ, even of m, in each individual is an indicator of the mood and physical status in which he or she is, either as a result of a sporadic situation, such as an episode of stress, a more severe disease or the age-related evidence for biological and physiological decline. We describe, in some detail, the most usual and clear methodology to calculate the phase space reconstruction because we have found that in its application to biological signals is not well understood at the physiological level, and its discriminant potential in the clinical setting no sufficiently exploit. In this sense, its effectiveness could be corroborated with the most modern state space reconstruction techniques that are less heuristic and with a more consolidated mathematical formalism 39 .