International Journal of Swarm Intelligence and Evolutionary Computation

Epidemiological Rogue Waves and Chaos-Induced Multifractal Self-Organized Criticality in COVID-19

Carlos Pedro Gonçalves^* Department of Management on Civil Aviation and Airports, Lusophone University of Humanities and Technologies, Lisbon, Portugal
^*Correspondence: Carlos Pedro Gonçalves, Department of Management on Civil Aviation and Airports, Lusophone University of Humanities and Technologies, Lisbon, Portugal, Email:

Received: 06-Apr-2024, Manuscript No. SIEC-24-25425; Editor assigned: 10-Apr-2024, Pre QC No. SIEC-24-25425 (PQ); Reviewed: 24-Apr-2024, QC No. SIEC-24-25425; Revised: 02-May-2024, Manuscript No. SIEC-24-25425 (R); Published: 10-May-2024, DOI: 10.35248/2090-4908.24.13.367

Introduction

The Severe Acute Respiratory Syndrome Coronavirus 2 (SARS- CoV-2) responsible for the Coronavirus Disease 2019 (COVID-19), was first identified in the city of Wuhan in the Hubei province in People’s Republic of China and led to the COVID-19 pandemic. Different research into the virus, by different authors applying chaos theory methods, uncovered complex nonlinear dynamics and markers of chaos in the main epidemiological series for different countries and regions [1-6].

In [5] we applied a combination of topological data analysis, machine learning and chaotic time series analysis methods and uncovered evidence of a form of stochastic chaos characterized by emergent low-dimensional noise resilient chaotic attractors underlying the number of new positive cases of COVID-19 per million and the number of new deaths from COVID-19 per million for the different world regions.

In the number of new positive cases of COVID-19 per million we found that the attractors were all three-dimensional, with positive Lyapunov exponents and a high predictability when using adaptive Artificial Intelligence (AI) based systems, employing topological information to predict the main target epidemiological series. The main topological properties of these attractors were analyzed, applying further topological data analysis methods. Similar methods were also successfully employed in [6] for the US and Canada daily hospital occupancies from COVID-19 data, also uncovering the existence of noisy low-dimensional chaotic attractors.

In [5], two patterns were identified in the epidemiological dynamics, one was the multiple wave pattern, which characterized mainly the Africa region, the other dynamics was the rogue wave-like pattern, which characterized the remaining world regions.

Rogue waves are a phenomenon in hydrodynamics that has been identified to occur in the sea as the formation of very large waves that correspond to waves that are at least twice as large as the significant wave height, rogue waves have been identified in large datasets with events including waves of 25 meters [7].

The occurrence of rogue wave-like dynamics in the epidemiological context as large exponentially rising epidemiological peaks in infections that then rapidly die out characterize a strongly turbulent epidemiological dynamics and a high-risk scenario for healthcare response strategies, since these peaks will lead to an accelerated flooding of hospitals and a saturation of hospital resources that may increase the death toll. It can occur in viruses that have a high transmissibility and a high number of asymptomatic cases.

The fact that this rogue wave dynamics is a dominant pattern in SARS-CoV-2’s epidemiological profile, as addressed in [5], may have roots in the virus’ characteristics, including the high transmissibility of new variants, with a possible large number of asymptomatic cases that may enhance the contagion dynamics and lead to rapidly accelerating outbreaks. China is one of the countries rogue wave patterns, with epidemiological peaks occurring mainly in urban regions, with large population numbers and high population density.

While one might consider the possibility that the rogue wave events in the COVID-19 numbers could be rare events, outside the normal epidemiological dynamics of the virus, this is, and however, an unlikely hypothesis given the nature of the virus and the predominance of this pattern, as reported in [5]. We find, in the present work, further evidence that this is not the case, by applying multifractal analysis and Smart Topological Data Analysis (STDA), combining topological data analysis with machine learning and chaos theory applied to the daily number of new confirmed cases of COVID-19 in China.

By applying the multifractal analysis methods, we uncover the presence of multifractal scaling in the epidemiological series, which, in complex systems, typically characterizes a form of self-organized criticality linked to turbulent signals, corresponding to Multifractal Self-Organized Criticality (MSOC), which has been reported to occur in financial data, characterizing financial turbulence [8-14].

We also uncover that this multifractal pattern is induced, in the case of China, by an emergent low-dimensional strange chaotic attractor in an emergent dynamics characterized by a three- dimensional phase space, we study this attractor’s topological structure, uncovering a link between the extreme rogue wave dynamics and a three-winged structure of the attractor, which provides for an important evidence of how the theories of chaos and Self-Organized Criticality (SOC) can be jointly applied to epidemiology, to study the link between emergent strange attractors in epidemiological dynamics as well as to effectively predict this dynamics, providing for both an explanatory ground for major risk patterns and for an opportunity to employ adaptive topological A.I. systems and attractor reconstruction methods for epidemiological prediction and risk management.

In section 2, we review the main concepts, materials and methods employed in the present work. In section 3, we present the results and in section 4 we provide for a final discussion on the implications of the present work for epidemiology, risk science and the relevance of chaos theory and the theory of SOC to these fields.

Materials and Methods

Self-organized criticality, also known as SOC, was proposed in [15,16] as an explanation of power law decay in power spectra of complex systems’ dynamics, growing as a theory within the complexity research field, to address the emergence of power law temporal and/or spatial correlations, leading to a scale invariance in time and/or in space as a result of complex systems’ self- organization.

The frequent occurrence of power law decay in power spectra, also called 1/f-noise, in different systems’ dynamics was the main problem addressed in [15], in this way, the hypothesis of SOC was initially proposed as an explanation for this phenomenon. The initial proposal of SOC, therefore, addressed specifically fractal scaling signatures, especially associated with power law scaling in spectral analysis.

Indeed, while scale invariance in power law decaying correlations occur at critical points of phase transitions, the difference in SOC is such that features like power law decay in spectral analysis, event size distributions or fractal laws in spatially extended systems result from a self-organization, rather than a careful adjustment of a parameter, in this way, the system is considered to self-organize to a critical dynamical regime, where scale invariance occurs. The theory therefore addresses a bottom-up rather than top-down occurrence of fractal laws in systems’ dynamics.

However, in some systems’ dynamics, the hypothesis of a single fractal scaling law holding with a single characteristic exponent, that is, a monofractal scaling does not hold, instead the emergent dynamics shows a spectrum of multiple scaling exponents, that is, the emergent scaling is multifractal rather than monofractal [17].

Some systems near critical points can exhibit multifractal scaling [18]. In this way, multifractal scaling can occur at critical points, which raises the possibility of SOC occurring with multifractal, rather than monofractal scaling, an issue that has been studied and identified in models of SOC and in dynamics of complex systems that range from financial returns’ series to earthquakes and even the Earth’s magnetotail [12,14,19-21].

A system’s self-organization to a multifractal dynamical regime corresponds to what can be called Multifractal Self-Organized Criticality (MSOC), which is an expansion of the theory of SOC to deal with the emergence of multifractal scaling in systems’ dynamics. In this case, the criticality is not characterized by a single scaling law but, instead, by a spectrum of scaling laws, that is, a multifractal spectrum.

Given a signal x(t) , the multifractal scaling is given by the moments’ order relation [17]:

log{E[|x(t)-x(t-s)|q]}~qh(q)log(s).........................................(1)

Keeping the lag fixed, the dependence upon the order is given by the product of the order by the generalized Hurst exponent h(q) which can be a nonlinear function of the moment order. Keeping the order fixed there is a power law scaling with the lag, thus, different orders lead to different scaling slopes, which explains the basis of the multifractal spectrum, that is, different power scaling laws (fractal laws) for different moment orders.

When the dependence of the Hurst exponents upon the order corresponds to a general nonlinear function, the logarithm of the moments of the absolute variations for each lag, have a complex nonlinear dependence upon the moments’ order, which can depend upon the structure of the multifractal spectrum. In this case, the multifractal scaling function addressed in [17] becomes relevant:

τ(q)=qh(q)-1........................................(2)

Depending upon the product of the moment order by the generalized Hurst exponent, the scaling function reveals the pattern for a possible moment divergence with the order or alternatively for a convergence. If we replace the multifractal scaling function in equation (1) we get the following relation [17]:

log{E[|x(t)-x(t-s)|^q]}~(1+τ(q))log(s).........................................(3)

Multifractal scaling is compatible with different power spectra, in this way, even if one has a white noise spectrum, multifractal scaling can still be present [8,10,17]. In this case, the emergence of multifractal scaling is a more complex problem than monofractal dynamics, in particular, as stressed in [17], spectral analysis can be misleading, in the sense that multifractal scaling can be present even in the case of white noise spectra, in this sense, the standard spectral analysis of a signal is insufficient to detect multifractal scaling, requiring the application of multifractal analysis methods for the estimation of the multifractal spectrum.

While multifractal signals can be built top-down with generators and multiplicative cascades, with these generators having played a key role in discussions around multifractal scaling and turbulence [17], the major issue, for both complexity research and risk science is that of self-organization of a system’s dynamics to a multifractal scaling regime, this is the already referred phenomenon of MSOC.

The presence of multifractal scaling can be analyzed through different methods, the method we will employ is Multifractal Detrended Fluctuation Analysis (MFDFA) with polynomial fitting, which is robust in the detection of monofractal versus multifractal scaling, and, in the case of turbulence processes, allows one to characterize the relation between risk and predictability. The method that we use is described in detail in [22], and it involves the estimation of a detrended fluctuation function that, for a fractal or multifractal signal, scales with the lag in accordance with a power law scaling rule that depends upon the generalized Hurst exponents, that is, for a multifractal process, the detrended fluctuation function F_q scales with the lag s as:

F_q~s^h(q) ............................(4)

Monofractal scaling can be identified when parallel lines in the log-log plot have the same slope, in this case, the generalized Hurst exponents coincide with each other or, at least, fluctuate in a small interval around the true Hurst exponent that characterizes the monofractal dynamics. Multifractal scaling is characterized by different slopes.

Besides the plot of the detrended fluctuation function for the different lags, we also plot the multifractal scaling function, the histogram for the generalized Hurst exponents and the graph of the function of the exponents for the different moments’ orders, besides these plots we also calculate the maximum and minimum of the generalized Hurst exponents’ distribution, allowing us to better characterize the profile of the multifractal scaling.

In the epidemiological context, MSOC resulting from the dynamics involving an epidemic or pandemic is a key point demanding further research. The epidemiological dynamics of a virus will depend upon biological factors including the incubation period, the first time for symptoms of a disease caused by the viral infection to appear, the number of asymptomatic cases, the viral mutation process and the mutation rates leading to the surfacing of new variants, the transmissibility and transmission vectors, the population it infects and immunity response of the human populations, with different factors, including human genetics, possibly affecting that response, the presence of biological reservoirs for the virus, other factors also include geographical and climate conditions, the country populations’ behaviors as well as public healthcare response policies and other measures including quarantine, possible treatment solutions, testing frequencies, the use of different measures, vaccination and so on.

Public healthcare response involves responses to predictions based on the epidemiological risk variables, which can lead to preemptive measures by healthcare authorities based on the predicted epidemiological dynamics; this introduces a feedback loop in the epidemiological process, with the adaptive responses impacting the epidemiological dynamics and vice-versa.

The above factors indicate that the epidemiological dynamics will tend to exhibit strong nonlinearities and feedback between top- down and bottom-up dynamics, especially in what regards healthcare planning and response as well as public health communication and the role of traditional and online media which may influence people’s behavior.

The above different factors will vary from country to country; this also includes the healthcare authorities and government approach to healthcare. For instance, in countries where the response is reactive rather than proactive in responding to rising healthcare pressure associated with viral infections, the dynamics will be different from the countries that implement proactive measures.

Considering the above points, the occurrence of multifractal turbulence in epidemiological dynamics, associated with MSOC, implies a complex process associated with this dynamics, in particular, in the case of rogue wave patterns, it involves possible sudden accelerations in the epidemiological processes, leading to extreme outbreak events that can put an accelerated pressure on healthcare resources.

Thus, in epidemiological turbulence, multifractal turbulence, associated with MSOC, leads to the possibility of sudden acceleration of outbreak processes leading to spikes of infections and to the possibility of epidemiological rogue waves of different sizes, that may occur at different times, with the largest rogue waves effectively appearing as very high spikes above the rest of the time series of confirmed cases.

From a healthcare standpoint this provides for a higher risk than the multiple wave patterns and the necessity for the identification of the source of MSOC. In [5], the rogue wave pattern was found to be predominant in SARS-CoV-2’s epidemiological dynamics and was also found to occur in tandem with chaotic markers, in this sense, one possible source of MSOC is a stochastic chaos dynamics, which is an open chaos dynamics.

Therefore, after testing for the presence of multifractal scaling and characterizing the scaling and memory properties of the COVID-19 epidemiological series, performing spectral analysis as a complementary analysis, we proceed to implement chaotic time series analysis methods, combined with STDA, in order to research the hypothesis of chaos-induced MSOC.

The first step is the time series phase space embedding for possible attractor reconstruction, where the reconstructed underlying phase point, for an embedding dimension d and a lag h, is given by:

p(t)=(x(t-(d-1)h),...,x(t-2h),x(t-h),x(t)).......................(5)

The main problem is one of finding the appropriate embedding dimension and embedding lag values.

The first parameter to be chosen is the embedding lag. As discussed in [5], in epidemiological contexts, the choice of the lag should be defined in terms of the incubation period and quarantine window, setting the period to the first day after a recommended or mandatory quarantine period allows for the embedding to account for possible quarantine effects.

In this way, as in [5], we use a 15-days period, which is the number of days between the start and the end of the World Health Organization’s (WHO) recommended a 14-day quarantine period from a person’s last exposure, leading to 15 days between the last exposure and an exit from quarantine.

As for the embedding dimension selection, we use STDA [5,6,23], employing an adaptive topological learner, that is, an artificial agent equipped with a topological machine learning unit, in this case we will use a k-nearest neighbors’ learning unit, the adaptive component involves employing a sliding learning window to predict the target epidemiological variable using the embedded phase point.

The use of a k-nearest neighbors learning unit is linked to the objective of the topological data analysis, in this case, we wish to find the optimal dimension, from a set of alternative dimensions, that leads to the highest exploitable topological information in the prediction of the target, in order to better capture the underlying attractor’s topology, to be able to apply topological analysis methods to characterize that attractor, linking the time series properties to that of the attractor’s structure. Since we will be using k-nearest neighbors’ graph analysis, with Euclidean metric, the machine learning unit should be a k-nearest neighbors’ unit [23].

The sliding learning window allows for an adaptation to local patterns in the attractor, including attractor epochs, making the agent adaptive to local topological regularities that can be exploited for the target prediction [5,6,23].

Thus, to select the embedding dimension we deploy the adaptive AI system, using a sliding window for relearning, in order to perform the single period prediction for the next period’s value of the target series, that is, the agent forms the following conditional expectation:

E[x(t+1)|p(t),w,k]=f(p(t),w,k,t).......................(6)

Where we denote by E[x(t+1)|p(t),w,k] the agent’s prediction for the next period’s target series value, this expectation is conditional on the phase point p(t) , on the learning window size w and on the number of nearest neighbors k.

The adaptive topological AI system is trained using the feature set of embedded points {p (t-w-1),...,p(t-2),p(t-1)} and target variable values {x(t-w),...,x(t-1),x(t)}. In this way, the AI learns to predict the next value of the target using the previous value of the embedded phase point. As stated, the relearning allows the AI to produce a local topological mapping, capturing the local topological regularities in the attractor and adapting to the attractor’s epochs, making the method especially robust in the case of turbulent and noisy series, which is effective when dealing with open systems and stochastic chaos as shown in [6,23], the method is also robust to bifurcations, as shown in [5].

The process of dimension selection that we will employ involves the simultaneous calibration of the number of nearest neighbors’ parameter and the embedding dimension. In this way, we use a grid-based calibration, recording the coefficient of determination score for the adaptive topological learner for different embedding dimensions in a range of alternative dimensions and different values of the number of nearest neighbors’ parameter in a range of alternative values.

In this way, we select the number of nearest neighbors and the embedding dimension that lead to the best performance in prediction of the target series. These are the parameters, within the searched range, which have the highest topological information exploitable by the adaptive topological agent, for which the phase point can be used in the prediction of the target series’ next period’s value.

For the best performer, besides the coefficient of determination score we report additional prediction performance metrics, including the correlation between the adaptive topological agent’s predictions and the target, the quotient of the Root Mean Squared Error (RMSE) by the data amplitude and the explained variance [23].

A high correlation, low RMSE/amplitude ratio, high explained variance and high coefficient of determination score indicate that there are topological patterns in the reconstructed trajectory that can be exploited for the prediction of the target series and that we can apply topological data analysis methods to better characterize the resulting reconstructed attractor.

The first method that we employ, for the reconstructed attractor is Eckmann et al,. method for the Lyapunov spectrum estimation [24]. If a convergence occurs in the spectrum with one or more positive Lyapunov exponents, this is indicative of chaos, a negative sum of exponents, with a positive largest Lyapunov exponent, is indicative of dissipative chaos, which is consistent with a chaotic attractor.

Using the estimated Lyapunov spectrum we also calculate the Kaplan-Yorke dimension [25,26], which is a measure of the fractal dimension of a chaotic attractor that uses the Lyapunov spectrum, the dimension calculation is based on the following formula:

Where, λ_j are the Lyapunov exponents arranged in decreasing order and n is the maximum number of exponents which may be added, with the exponents arranged in decreasing order, before the sum becomes negative.

Besides the Kaplan-Yorke dimension we also calculate the Lyapunov time, which is the inverse of a dynamical system’s largest Lyapunov exponent and sets the limits of predictability.

As a complementary analysis to the Lyapunov time we calculate the predictability of the target series several days ahead, using the adaptive topological agent and the embedded trajectory for prediction, thus, instead of the one-day ahead prediction of equation (6), we use the same sliding window adaptive topological learner scheme for an l-days ahead prediction, using as a predictor the phase point at time t, leading to the prediction values.

We increase the number of days ahead for prediction and, in each case, extract the corresponding coefficient of determination score and report the first time that the score drops below a threshold, we call this the threshold foresight time. In this case, we report the 60% and the 50% foresight times. The foresight time corresponds to the time it takes for the exploitable topological patterns for long- term prediction by the adaptive topological learner to drop below a threshold. We analyze the relation between the foresight time marks and the estimated Lyapunov time.

After analyzing the Lyapunov spectrum, the Kaplan-Yorke dimension, the Lyapunov time and the foresight times, we address the topological features of the attractor and specific computational properties associated with symbolic dynamics’ analysis by performing three analyses: The k-nearest neighbors’ graph analysis, persistent homology analysis and permutation entropy analysis calculated from the permutation graphs.

The k-nearest neighbors’ graph, which links each phase point in the reconstructed attractor to each of its k nearest neighbors is an important topological analysis method for the characterization of the attractor’s topological structure, we begin by plotting the graph and the degree distribution for the number of nearest neighbors for which we obtained the best prediction performance, we, then, vary this number and, for each resulting graph, we calculate the following metrics:

• The percentage of completeness: Defined as the quotient of the graph’s number of edges over the total number of edges that would make the graph complete multiplied by 100, this metric provides for a measure of the degree to which the graph is close to a complete graph, expressed as a percentage;

• The graph’s normalized degree entropy: Defined as the Shannon entropy of the degree distribution divided by the maximum entropy, this entropy value is between zero and one, and the closer it is to one the closer the graph’s degree distribution is to a maximum entropy distribution;

• The graph’s Kolmogorov-Sinai (K-S) entropy: The entropy, for an unweighted graph, which in our case coincides with the logarithm of the dominant eigenvalue of the graph’s transition matrix, using the base 2 for the logarithm, it provides for an information measure, in this case, considering a random walker that follows the k-nearest neighbor’s adjacency, the K-S entropy graph provides for a measure of the rate at which information is generated by the graph’s adjacencies, providing, in this way, for a measure of the topological complexity of the attractor associated with the k-nearest neighbors’ structure.

We calculate each of these metrics for increasing values of number of nearest neighbors, studying how these metrics change with this number.

After performing this analysis, we proceed to plot the attractor in three-dimensional Euclidean space in order to have an overall picture of the three-dimensional projection of the attractor, in this case, as we will see, the attractor for China is three-dimensional and its overall geometric structure will give us significant information on the source of the rogue wave epidemiological dynamics.

Besides this plot, we also apply persistent homology analysis, which provides for a picture of persistent structures in the attractor, by which different homology classes are counted, for different homology dimensions, this analysis was applied in [5] to regional COVID-19 new cases per million and new deaths per million series.

Homology classes are classified by the dimensionality of the boundary of the structure, in this way, a zero-homology class corresponds to connected components in the attractor, therefore, with a zero-dimensional boundary, while a one-homology class corresponds to components that have a one-dimensional boundary, and so on. In persistent homology analysis the number of structures in each simplicial complex in a Vietoris-Rips filtration are counted for each homology class, the birth and death of homology classes as the radius is increased is also counted [5].

The lifetime or persistence of a class c , at a given homology dimension, corresponds to the difference between the death and birth times, with the death happening after birth, in this way, denoting the time of death by n_D and the birth time by n_B , the persistence of a homology class can be calculated as:

Structures that live through the full filtration have n_D = ∞ and, therefore, pers (c) = ∞ . In terms of statistical analysis of the homology classes, we calculate the following metrics [5]:

• The number of classes with pers(c) < ∞ : Calculated for each dimension, this metric provides for a statistical distribution of the number of structures for each homology dimension with lifetimes shorter than infinity;

• The number of classes with pers(c) = ∞ : This metric allows us to identify the homology dimensions that have structures that persist throughout the whole filtration, allowing us to identify large-scale structures;

• The maximum persistence: Calculated for each dimension, this metric allows us to characterize which dimensions have the strongest persistent structures;

• The mean persistence: Also calculated for each dimension, this metric allows us to characterize each homology dimension in terms of its mean persistence.

The last analysis method that we employ is a symbolic dynamics analysis using the ordinal partition graph, which is used to calculate the permutation entropy, this analysis was also employed in [5], but is expanded here to deal with patterns associated with different embedding dimensions.

Indeed, given the permutations Π_d{ π_i: i=1,2,......,d!} of the embedding dimensions’ set {1,2,.....,d} with i=1,2,...,d! one defines a permutation map:

That is, the permutation produces a non-decreasing reordering of the values in the d tuple of real numbers.

The permutation map defined above is calculated for each phase point in the embedded trajectory, leading to a sequence of permutations, this sequence leads to an ordinal partition graph, such that two permutations are linked if there is a transition between them. From the ordinal partition graph, the permutation entropy can be calculated from the estimated probability distribution over the permutations as follows [5]:

Besides the permutation entropy, we also calculate the permutation graph’s K-S entropy, degree entropy and percentage of completeness, as we did for the k-nearest neighbors’ graphs.

For increasing embedding dimensions we obtain the corresponding permutation graphs and calculate the above metrics, leading to a picture of how these metrics change with the embedding dimension, allowing us to analyze the correlations between each of these metrics with the coefficients of determination obtained for the best performer adaptive topological learner selected in the calibration process, that is, using the number of nearest neighbors that leads to the highest coefficient of determination in the calibration process for the optimal embedding dimension, we study the relation between the prediction performance of this agent with varying dimensions and the ordinal partition graph’s metrics, thus providing for a better picture of the relation between the computational properties associated with the ordinal partition graphs and the prediction performance with the embedding dimension.

The above analysis methods are applied to the number new confirmed cases of COVID-19 in China, from 2020-01-03 to 2023- 10-27. Allowing us to characterize the underlying emergent attractor and uncover the origin of the rogue wave dynamics.

Results and Disussion

In Figure 1, we show the time series chart for the daily new confirmed cases of COVID-19 in China in the period from 2020-01-03 to 2023-10-27, obtained from the “Our World in Data” database. As can be seen from the graph, it almost seems as if there is a flat line and then some peaks, near the end with a large wave of confirmed cases, this is actually an illusion of scale, indeed, the region that appears as a flatline is actually comprised of fluctuations in the number of new cases, and large jumps as illustrated in Figure 2.

Figure 1: Daily new confirmed cases of COVID-19 in China in the period from 2020-01-03 to 2023-10-27.

Figure 2: Figure 1’s series from: 2020-01-03 to 2021-04-16 (left), 2021-04-17 to 2022-07-30 (middle), 2022-07-31 to 2023-10-27 (right).

The two figures exhibit a pattern of turbulence in the series with extreme jumps corresponding to the epidemiological rogue wave dynamics of different sizes occurring at different steps, however, this dynamics follows the same rules of structure as the whole process, as we now show applying fractal and multifractal analysis tools.

In Figure 3, we show the power spectrum plotted in log-log scale (left) and with the fitted line in the power law decaying region (right). As can be seen in the figure, the power-law scaling occurs in the high frequency region, the estimated slope is of -3.9562 with an R² of 99.69%, which is evidence favorable to a strongly persistent process in the black noise range, that is, a power law noise with decay exponent larger than 2, a feature that was also identified in the regional data in [5].

Figure 3: Power spectrum plotted in log-log scale (left) and with fitted line (right).

Now, while the data indicates long-memory persistent process, which is consistent with SOC, we need to analyze the multiscaling properties in order to evaluate the hypothesis of MSOC. In this case, we apply Multifractal Detrended Fluctuation Analysis (MFDFA) to estimate the multifractal spectrum and the extract the main multifractal metrics (Figure 4).

Figure 4: Multifractal detrended fluctuation analysis, showing the: fluctuation function (top left), generalized Hurst exponents (top right), scaling function (bottom left) and generalized Hurst exponents’ histogram (bottom right), the moments’ orders range from 0 to 40 in 200 steps and the lags range from 1.5 to 2.5 in 500 steps.

As shown in Figure 4, we find a few important points, for moments’ orders ranging from 0 to 40, there is a wide range of slopes which is indicative of multifractal scaling (Figure 4, top left), also, the estimated exponents exhibit a nonlinear dependence upon the order, with an accelerated decay to exponent values just below 0.5 (Figure 4, top right), which leads to a divergence in the moments as can also be seen in the multifractal scaling function (Figure 4, bottom left) which has an initial nonlinear rise and then rises linearly with the order, which is not consistent with a moment convergence with the order.

The histogram reflects the convergence of the exponents to the values near 0.5, indeed, the dominant exponents, in terms of frequencies, are near 0.5 (Figure 4, bottom right), with the maximum estimated exponent being around 2.5196 and the minimum estimated exponent being around 0.4987.

The relation of the scaling with regards to the order and lag is relevant, in this case, for the high moments’ orders there is a convergence to a monofractal scaling with a Hurst exponent below but near 0.5 with respect to the lag, however, for the low moments’ orders, the Hurst exponents are higher than 1, which is significant.

The fact that the Hurst exponents for the lower moments’ orders are greater than 1, implies that there is a divergence of the lower moments with the lag due to higher than 1 exponent, which is consistent with the high turbulence process associated with the dynamics of the daily new confirmed cases of COVID-19.

The overall results of the multifractal analysis indicate that the dynamics is multifractal, rather than monofractal and the spectrum already shows a relation to the turbulence features. This is evidence favorable to the hypothesis of MSOC. Also of notice, the multifractality occurrs with a power spectrum that is not white, that is, we have a black noise spectrum for the signal, which indicates strong persistence, along with a multifractal scaling structure.

In terms of epidemiological dynamics, MSOC associated with the turbulence pattern, means that the epidemiological dynamics associated with the rogue wave pattern in China is indicative of a self-organization to a critical multifractal turbulent regime where the large jumps are part of the same process that leads to the smaller jumps. The question that needs to be raised is whether the source of the MSOC is an underlying chaotic dynamics, emergent from the epidemiological dynamics, that is, whether we are dealing with chaos-induced MSOC.

As discussed in the previous section, using the 15 day period for delay embedding, we search for the embedding dimension with the highest exploitable topological information in the prediction of the daily number of new confirmed cases of COVID-19 in China, using the adaptive topological learner.

In Table 1, we present the results for the embedding dimensions ranging from 2 to 10 and the nearest neighbors’ parameter k ranging from 2 to 6. Three patterns stand out in the results; the first is that the prediction performance is very high for all embedding’s, with coefficients of determination scores all greater than 79%, which indicates a strong topological exploitable pattern for the prediction of the target series. The second pattern is that for all values of k, the best performance is obtained for a three-dimensional embedding, which is strong evidence of an emergent attractor with three degrees of freedom underlying the epidemiological dynamics. A point that is particularly relevant since in [5] we also found three- dimensional attractors for all the world regions that had a stable attractor structure, which are all world regions with the exception of Oceania which underwent a bifurcation.

Table 1: R² score for the adaptive topological learner in predicting the daily number of new cases of COVID-19 in China with a k-nearest neighbors’ learning unit, a 15-days sliding learning window, using the previous day’s phase point as predictor, with embedding dimensions ranging from 2 to 10, k ranging from 2 to 6 and weights based on Euclidean distance, with ball tree algorithm used for the learning process.

The third pattern that we find is that the performance decreases with the increase in the number of nearest neighbors, a pattern that is different from the study by continent in which the performance increased with the number of nearest neighbors [5].

The best performance of the group is obtained for k=2 and a three-dimensional embedding, which leads to a coefficient of determination score of 87.567%.

In Table 2, we show the additional prediction metrics for k=2 and d=3. We find, in this case, that the correlation between the adaptive agent’s predictions and the target series is positive and high, with a value of 0.9393, which reinforces the results on a good prediction performance using the k-nearest neighbors algorithm, that is, the k-nearest neighbors for the reconstructured attractor contains exploitable information that can be used for predicting the target series with a high prediction performance.

Table 2: Prediction performance metrics for table 1’s best performing adaptive topological learner of the previous table, with embedding dimension equal to 3.

The RMSE is only 2.615% of the total amplitude of the data, which is a low relative error; the explained variance and the coefficient of determination score are both slightly higher than 87.5%. In Figure 5, we show the time series chart with AI’s predictions and the daily number of confirmed cases of COVID 19 in China.

Figure 5: Daily new confirmed cases of COVID-19 in China in the period from 2020-01-03 to 2023-10-27 in blue and prediction by the best performing adaptive topological learner in orange.

The Figure illustrates visually the results that the statistics provide, that is, a high prediction performance for the topological learner using the three-dimensional embedding and the 15-day period lag. Thus, the results obtained from the topological learner lead to strong evidence in support of the possibility of a noise-resilient three-dimensional attractor, underlying the epidemiological dynamics.

The next step in the analysis is the estimation of the Lyapunov spectrum, in this case, we find a convergence of the maximum Lyapunov exponent to a positive value, up to a four-decimal place approximation, of 0.0322, as shown in Figure 6, with the other two exponents being negative and, respectively, up to a four-decimal place approximation equal to -0.0036 and -0.0531, which leads to a negative spectrum sum of -0.0244, the convergence of the spectrum, with positive largest exponent and negative sum are consistent with dissipative chaos. In this case, the Lyapunov time is, up to a four- decimal place approximation, of 31.0181, which is between 31 and 32 days, the Kaplan-Yorke dimension is, in turn, 2.5393 which is between 2 and 3 dimensions.

Figure 6: Lyapunov spectrum estimation for the reconstructed three- dimensional attractor.

The high predictability with best results showing consistently a three-dimensional embedding as the best embedding in the grid search, the convergence of the estimated Lyapunov spectrum to a positive largest Lyapunov exponent, the negative sum of the spectrum and the Kaplan-Yorke dimension are all consistent with a noise-resilient low dimensional strange chaotic attractor underlying the target series’ dynamics. This evidence and the previously found markers of MSOC are favorable to the hypothesis of chaos-induced MSOC.

The largest positive Lyapunov exponent along with the strongly persistent process with a black noise signal identified in the power spectrum analysis of the series and a multifractal scaling of this same series leads to an added complexity in regards to the predictability horizon.

Indeed, the positive Lyapunov exponent leads to an exponential divergence and restricts the predictability, in this case to a Lyapunov time that is between 31 and 32 days, in this way, the loss in predictability associated with the sensitive dependence upon initial conditions is expected to be identified at 32 days, however, the long memory of the dynamics and multifractal scaling lead to topological predictability markers that may increase predictability.

In order to address the pattern of predictability we now calculate the coefficient of determination for the adaptive topological learner, with the three-dimensional reconstructed attractor in several periods’ ahead predictions, and calculate the foresight time for 60% and 50% prediction, this is shown in Figure 7.

Figure 7: Prediction performance of the best performer in table 1, for t-periods ahead prediction, with t varying from 1 period (1 day) to 50 periods (50 days) using the three-dimensional embedding, 15 days’ lag and the R2 score as prediction performance measure. Note: () R²; () 60% line; ()50% line.

We find that the predictability decreases with the prediction horizon, when we use the best performing topological learner to predict several periods ahead, rather than one period ahead, with the prediction performance measured by the coefficient of determination first dropping below 60% from 24 to 25 days ahead, making the 60% foresight time 25.

The 50% foresight time, in turn, is 32. with the coefficient of determination first drop below 50% occurring from 31 to 32 days ahead prediction, which coincides with the expected foresight from the Lyapunov time, which is also between 31 and 32 days, a result and reinforces the strong evidence of a low-noise three-dimensional chaotic attractor.

There is, however, a plateau with the increasing of the window in which the predictability seams to stabilize fluctuating in a band between 60% and 50% for the coefficient of determination score, which may be linked to the topological features associated with the long memory, the fluctuations in the narrow band last up until a prediction horizon of 46 days ahead, after this period at 47 days the performance drops and then accelerates with the coefficient of determination decaying significantly only after the 46 period horizon, which is around 1.4830 of the Lyapunov time.

Another relevant point is that the predictability of the topological adaptive agent does not decrease exponentially until the foresight time of 32 days ahead but, instead, decreases linearly, this may be linked to the attractor’s structure, namely, the topological information in the k-nearest neighbors and the long memory persistence of the series associated with the MSOC with a black noise power spectrum which leads to the strong persistence.

The results so far are consistent with a form of chaos-induced MSOC, with an underlying low-noise three-dimensional chaotic attractor, with a strong topological structure that allows for a high- predictability of the target series and with the predictability limits linked to the Lyapunov time and the predictability pattern related to the interplay between the MSOC with black noise spectrum (black MSOC) and the positive largest Lyapunov exponent.

Having identified the presence of an attractor we now need to apply topological analysis methods in order to better characterize its structure and its relation to the predictability of the target series. We begin by addressing the topological patterns based upon the number of nearest neighbors. In Figure 8, we show the k-nearest neighbors’ graph for 2 nearest neighbors and the corresponding degree distribution plotted in log-log scale. For this graph, the percentage of completeness is low, around 0.201%, the degree entropy is also low 0.1534, indicating a far from maximum entropy degree distribution, the K-S entropy is, in turn, equal to 3.0543 bits. Regarding the degree distribution itself, the graph does not show a scale-free structure but, instead a faster than power law decay.

Figure 8: k-nearest neighbors’ graph for 2 nearest neighbors (left) and respective degree distribution (right) plotted in log-log scale, for the reconstructed attractor.

This distribution profile holds for higher number of nearest neighbors, however, there are power law scaling relations between the graph’s number of nearest neighbors and the graph’s main metrics, and the predictability as shown in Figure 9. First, considering the percentage of completeness of the graph, defined, as reviewed in the previous section, in terms of the relation between the number of edges in the graph and the number of edges that would make the graph complete, we find that the percentage of completeness grows as a power law of k, the same holds for the degree entropy and the K-S entropy, while the coefficient of determination score for the adaptive topological agent decays as a power law of k as shown, for k varying between 2 and 12.

Figure 9: Power law scaling of the percentage of completeness, degree entropy, K-S entropy and R² score for the adaptive topological agent in the one period ahead prediction, with k varying from 2 to 12 in steps of 1.

The above results provide for a link between the topological structure of the attractor and the predictability, namely, the main graph’s metrics increase as a power law of the number of nearest neighbors with the predictability decreasing as a power law. In this sense, the power law rise in the graph’s metrics, in particular the entropy measures that provide for a definition of graph complexity, explains the power law decay in the predictability with the rise in the parameter k.

In this sense, as k is increased, the exploitable topological information decreases due to a rise in the complexity of the neighborhood connectivity structure addressed in terms of the k-nearest neighbors.

So far, the evidence is favorable to a chaos-induced MSOC, associated with the daily number of new cases from COVID-19 in China, with an emergent noise-resilient low-dimensional attractor underlying the turbulent dynamics. The k-nearest neighbors’ analysis combined with the chaotic time series analysis methods effectively allowed us to link the Lyapunov time and the main graph metrics (the percentage of completeness, the degree entropy and the K-S entropy) associated with the k-nearest neighbors’ graph to the predictability of the series by the adaptive topological learner.

Having performed the k-nearest neighbors’ topological analysis, we now turn to the persistent homology analysis. In Figure 10 (left), the reconstructed attractor is shown. As can be seen the attractor exhibits a three-winged structure, this structure accounts for the rogue wave phenomenon with the tips of the wings leading to the rogue waves. In Figure 10 (right), we show the persistence diagram for homology dimensions 0 and 1. From the diagram we find that all homology dimensions 0 structures are born at the beginning of the filtration exhibit the strongest persistence, with an infinity class, for homology dimension 1, which represent cycles. We find that most classes are short-lived with just one class lasting longer (Figure 10).

Figure 10: Reconstructed attractor on the left and persistence diagram on the right. Note: () ∞; () H₀ ; ()H₂.

The results from Figure 10 are reinforced in Table 3, which shows the persistence statistics.

Table 3: Persistent homology statistics for the reconstructed attractor.

In terms of persistence, as shown in Table 3, we find that the homology dimension 0 is predominant with one infinity class and 1344 non-infinity classes, all born at 0, the maximum persistence of the non-infinity classes is 1008368.06 and the mean persistence is 25166.92, while homology dimension 1, has only 375 classes, higher maximum persistence (1876963.88) but lower mean persistence 5856.45.

The persistence values contrast significantly with those found for the Asian continent’s chaotic attractor in [5], which was also three-dimensional but had a lower number of homology dimension 0 and homology dimension 1 classes (899 in for homology dimension 0 and 257 for homology dimension 1) with lower persistence statistics, which indicates the presence of a greater number of topological structures for China’s attractor with stronger periodicities.

Considering now the ordinal partition analysis, obtained for the reconstructed attractor, Figure 11 shows the ordinal partition graph, the graph is complete, which means that transitions between all possible alternative permutations occur.

Figure 11: Ordinal partition graph for the three-dimensional reconstructed attractor.

The permutation entropy is up to a four decimal places’ approximation equal to 2.5023 bits, the degree distribution entropy is 0, and because the graph is complete the K-S entropy of the graph is 2.3219.

The permutation analysis becomes more important in terms of characterization of the dynamics by considering embedding’s with increasing embedding dimension. In Figure 12, we show the different entropy measures, permutation entropy, K-S entropy, degree entropy for the ordinal partition graph, and the percentage of completeness as a function of the embedding dimension for embedding dimensions varying from 2 up to 8.

Figure 12: Permutation entropy (top left), K-S entropy (top right), degree entropy (bottom left) and percentage of completeness (bottom right) plotted for embedding dimensions varying between 2 and 6 for the ordinal partition graphs.

Considering the results from Figure 12, we find several relevant points. First, the percentage of completeness shows that for 2 and 3 dimensions the permutation graphs are complete, that is, the embedded trajectory leads to a computational dynamic where all of the alternative transitions between permutations occur, this also explains why for 2 and 3 dimensions the degree entropy is 0.

The ordinal partition graph’s K-S entropy, however, increases when we go from 2 to 3 dimensions and increases again when we go to a four-dimensional embedding. The degree entropy also increases when we proceed from a three-dimensional to a four- dimensional embedding, in this case, not all of the transitions between permutations are allowed, indeed, in 4 dimensions, the percentage of completeness drops, in this case, to 69.203%. The percentage of completeness continues to drop with the embedding dimension until it reaches almost zero, for 7 and 8 dimensions. This is so because the dynamics is such that number of allowed transitions grows slower than the number of possible transitions between permutations.

In terms of ordinal partition graph structure the degree entropy is always lower than the maximum entropy, which means that we do not have an equiprobable distribution of the degree distribution, furthermore the largest entropy is achieved for 5 dimensions, the second largest being achieved for 4 dimensions. With increasing embedding dimension beyond 5 there is a rapid decrease in degree entropy, which makes the degree distribution far from the maximum entropy. The permutation entropy rises linearly with the embedding dimension; however this rise slows down in the transition from 7 to 8 dimensions.

The most relevant metric, in this case, is the K-S entropy, since not only does it provide for a measure of the rate at which information is generated for the permutations’ transition graph, its pattern shows a specific relevance in the relation with the exploitable topological information, namely, it is positively and strongly correlated with the coefficient of determination score obtained for the adaptive topological agent with k equal to 2, as shown in Table 4, we find that the correlation is, in this case, 0.8706, which is positive and strong, that is, the higher the K-S entropy for the permutation graph tends to be, the higher tends to be the performance of the adaptive topological learner. The remaining entropy measures (permutation entropy and degree entropy) are also positively correlated with the R² score, however, the correlation is weaker, the percentage of completeness, on the other hand, has a negative but also weak correlation with the prediction performance.

Table 4: Correlations between each metric for the permutation graphs obtained for the different embeddings and the corresponding R² scores obtained for the adaptive topological learner with 2-nearest neighbors’ learning values provided in Table 1.

In terms of symbolic computational analysis, considering the above results, we find that the graphs for the permutation transitions obtained from the embedded trajectories are such that the K-S entropy of these graphs shows strong links to the exploitable topological information in the embedded trajectory for predicting the target, with varying embedding dimensions, thus, the K-S entropy of the permutation graphs stands out in the sense that the higher its value, the higher tends to be the predictability of the series when using the embedded trajectory as predictor feature space by the adaptive topological learner.

This pattern is linked to the fact that, as the embedding dimension is increased, the lower tends to be the K-S entropy and also the lower tends to be the predictability of the series, in terms of exploitable topological information. This in turn reinforces the evidence for low-dimensional noisy chaotic dynamics.

Conclusion

The epidemiological series for COVID-19 has a strongly turbulent process with rogue waves as a predominant pattern, in the case of China, the lockdown policies actually led to a curbing of the wave heights, however, rogue waves occurred at different heights, strongly linked to urban centers’ transmission and the rise of new variants, which, coupled to the virus’ epidemiological profile, characterized by high transmissibility and large number of asymptomatic cases, leads to a fast propagation dynamics in human populations.

We found that this rogue wave turbulent dynamics has markers of Multifractal Self-Organized Criticality (MSOC), with a black noise spectrum linked to a strong persistence.

The application of state-of-the-art chaos theory empirical methods employing Smart Topological Data Analysis (STDA) uncovered evidence that the MSOC in the daily number of new confirmed cases of COVID-19 in China is linked to an emergent three- dimensional noise resilient chaotic attractor with a three-winged structure, with the tips of the wings in phase space explaining the rogue wave phenomenon. The signal is strongly predictable from the reconstructed attractor by an adaptive topological learner using a k-nearest neighbors’ topological unit, which can be used by healthcare authorities for response planning.

In several days ahead prediction we found a link between the Lyapunov time and the time it takes for the coefficient of determination score of the best performing adaptive topological agent to drop below 50% (foresight time), which further reinforces the findings for an emergent chaotic attractor.

Further topological data analysis on the attractor’s structure found a link between the topological structure of the attractor and the predictability pattern of the new cases of COVID-19 by the adaptive topological learner, with the degree entropy and the K-S entropy of the k-nearest neighbors’ graph strongly linked to the prediction performance of the topological learner, with varying values of the number of nearest neighbors’ parameter k.

The persistent homology analysis showed that the attractor has a higher persistence than the results obtained in a previous study for the Asian continent [5], which also had a three-dimensional emergent attractor.

The ordinal partition analysis methods uncovered a link between the computational dynamics associated with the phase space embedding, for different embedding dimensions, and the predictability of the target series with the K-S entropy of the ordinal partition graph being the main metric positively and strongly correlated with the predictability of the series, accounting for the drop in predictability as the embedding dimension rises beyond the three-dimensional structure.

The overall strong evidence for chaos-induced MSOC and the strong predictability of the target series, employing an adaptive topological learner, constitutes a relevant finding for both the strategic planning for dealing with COVID-19 in countries and regions with rogue wave patterns and for epidemiological studies, since it provides an example where there is a robust evidence for an emergent chaotic attractor in epidemiological dynamics, in the context of a pandemic.

From an epidemiological and risk science standpoint, the evidence indicates the importance of further research into emergent chaos-induced self-organized criticality in the COVID-19 as a case study for possible future pandemics, also illustrating the need for the application of latest methods from chaos theory combining topological data analysis and machine learning, in order to understand the emergence of low-dimensional attractors underlying epidemiological turbulence markers with multifractal scaling of target epidemiological variables and the link between the epidemiological dynamics and the attractor’s geometric and topological properties.

The results also reinforce the effectiveness of the employment of adaptive A.I. systems with machine learning units and sliding windows for target epidemiological series’ prediction and healthcare response planning.

References

1. Mangiarotti S, Peyre M, Zhang Y, Huc M, Roger F, Kerr Y. Chaos theory applied to the outbreak of COVID-19: An ancillary approach to decision making in pandemic context. Epidemiol Infect. 2020;148:e95.

   [Crossref] [Google Scholar] [Pubmed]

2. Postavaru O, Anton SR, Toma A. COVID-19 pandemic and chaos theory. Math Comput Simul. 2021;181:138-149.

   [Crossref] [Pubmed]

3.  Sivakumar B, Deepthi B. Complexity of COVID-19 dynamics. Entropy. 2022;24(1):50.

   [Crossref] [Google Scholar]

4. Borah M, Gayan A, Sharma JS, Chen Y, Wei Z, Pham VT. Is fractional-order chaos theory the new tool to model chaotic pandemics as Covid-19? Nonlinear Dyn. 2022;109(2):1187-1215.

   [Crossref] [Google Scholar] [Pubmed]

5. Gonçalves CP. Low dimensional chaotic attractors in SARS-CoV-2’s regional epidemiological data. Int J Swarm Evol Comput. 2022;11(9): 1000271.
6. Gonçalves CP (2023). Low dimensional chaotic attractors in daily hospital occupancy from covid-19 in the usa and canada. Int J Swarm Evol Comput, 2023;12(1):1000291.
7. Nikolkina I, Didenkulova I. Rogue waves in 2006–2010. Nat Hazards Earth Syst Sci. 2011;11(11):2913-2924.

   [Google Scholar]

8. Mandelbrot BB, Fisher A, Calvet L. A multifractal model of asset returns. Cowles Foundation Discussion Papers. 1997:1164.

   [Google Scholar]

9. Mandelbrot BB. The variation of certain speculative prices. In Fractals and scaling in finance. Springer. 1997:371-418.

   [Google Scholar]

10. Mandelbrot BB, Hudson RL. The (mis) behaviour of markets: A fractal view of risk, ruin and reward. Basic Books. 2004.

    [Google Scholar]

11. Voit J. The statistical mechanics of financial markets. Springer. 2003.

    [Google Scholar]

12. Gonçalves CP. Financial turbulence, business cycles and intrinsic time in an artificial economy. Algor Fin. 2011;1(2):141-156.

    [Google Scholar]

13. Gonçalves CP. Quantum financial economics-Risk and returns. J Syst Sci Complex. 2013;26(2):187-200.

    [Google Scholar]

14. Gonçalves CP. Coupled stochastic chaos and multifractal turbulence in an artificial financial market. Int J Swarm Evol Comput, 2022 11(7): 1000261.
15. Bak P, Tang C, Wiesenfeld K. Self-organized criticality: An explanation of 1/f noise. Phys Rev Lett. 1987;59(4):381-384.

    [Crossref] [Google Scholar] [Pubmed]

16. Bak P, Paczuski M. Complexity, contingency, and criticality. Proc Natl Acad Sci.1995;92(15):6680-6696.

    [Crossref] [Google Scholar] [Pubmed]

17. Mandelbrot BB. Fractals and scaling in finance-discontinuity, concentration, risk. Springer. 1997.

    [Google Scholar]

18. Zhao L, Li W, Yang C, Han J, Su Z, Zou Y. Multifractality and network analysis of phase transition. PLoS One. 2017;12(1): e0170467.

    [Crossref] [Google Scholar] [Pubmed]

19. Pinho R, Andrade RFS, Nogueira Jr. E. Multifractal aspects for a self-organized critical model. Int J Mod Phys B. 2003;17:3065-3073.

    [Crossref] [Google Scholar]

20. Rahimiâ??Majd M, Shirzad T, Najafi MN. A selfâ??organized critical model and multifractal analysis for earthquakes in Central Alborz, Iran. Sci Rep. 2022;12:8364.

    [Crossref] [Google Scholar] [Pubmed]

21. Chang T. Self-organized criticality, multi-fractal spectra, sporadic localized reconnections and intermittent turbulence in the magnetotail. Phys Plasmas. 1999;6(11): 4137-4145.

    [Google Scholar]

22. Gorjão LR, Hassana G, Kurths J, Witthauta D. MFDFA: Efficient multifractal detrended fluctuation analysis in python. Comp Phys Commun. 2022;273:108254.

    [Crossref] [Google Scholar]

23. Gonçalves CP. Chaos-induced self-organized criticality in stock market volatility-an application of smart topological data analysis. Int J Swarm Evol Comput. 2023;12(5):331.
24. Eckmann JP, Kamphorst SO, Ruelle D, Ciliberto S. Liapunov exponents from time series. Phys Rev A Gen Phys. 1986;34(6):4971- 4976.

    [Crossref] [Google Scholar] [Pubmed]

25. Frederickson P, Kaplan J, Yorke E, Yorke J. The lyapunov dimension of strange attractors. J Diff Eqs. 1983;49(2):185–207.

    [Crossref] [Google Scholar]

26. Kaplan J, Yorke J. Chaotic behavior of multidimensional difference equations. In Peitgen, HO; Walther, HO (eds.). Functional Differential Equations and the Approximation of Fixed Points. Lecture Notes in Mathematics. Springer. 2006;730:204-227.

    [Google Scholar]