Introduction

It is possible to study any physiological process through an intricate network of biochemical reaction mechanisms from the systems biology paradigm^1,2,3,4, such mechanisms are responsible for regulating a wide variety of processes of utmost importance for life, through complex positive and negative feedback systems^5,6. A prominent example of negative feedback is the thyroid hormone regulatory mechanism⁷ carried out by either thyroid cells or the leptin-insulin axis⁸, which is strongly related to metabolic processes. Discrepancies in these feedback mechanisms can lead to multiple pathologies at the systemic level⁷. With the development of experimental tools for studying such complex systems, the need to understand the underlying dynamics of these regulatory mechanisms also arose, so biochemists undertook the task of studying the chemistry of these mechanisms to determine the kinetic parameters of the reactions participating in these processes. Meanwhile, biophysicists had to translate these biochemical processes to mathematical models using tools based on Dynamical Systems Theory (DST) and the discoveries made by biochemists^9,10,11.

B. Belousov was a Russian biophysicist who pioneered in the study of complex regulatory mechanisms present in biochemical processes at an experimental level. In 1950 he was given the task of proposing a reaction mechanism analogous to the Krebs cycle to study feedback processes^12,13,14. Belousov conceived a chemical mixture made essentially of citric acid, bromate ions, and cerium ions in an acid medium under constant agitation¹³. During the reaction, Belousov observed that the mixture changed from being transparent to a yellow hue and vice versa. The phenomenon occurred time after time, being reminiscent of the Krebs cycle feedback processes. Nevertheless, his work was never officially published because reviewers concluded that the process was caused by mixture impurities, and that the phenomenon violated the natural laws of thermodynamics^13,14. Years later, Prigogine defined the basis of thermodynamics of irreversible processes ^12,13,14.

Later on, another Russian biophysicist, A. Zhabotisky resumed Belousov’s work. He replaced citric acid with malonic acid and cerium ions with iron ions to visualize chemical species concentration-oscillations obtaining a solution that shifted from red to blue and vice versa¹². This chemical mechanism was thus named the Belousov-Zhabotinsky (BZ) reaction^12,13,14. The work of Belousov and Zhabotinsky set the basis for the study of oscillatory biochemical processes¹⁴.

In 1972, Field, Köros and Noyes proposed the first mathematical model in differential equations that described the underlying dynamics of the BZ reaction (FKN Model) and laid the groundwork for mathematical modeling of chemical mechanisms with oscillatory behaviors^15,16. This opened a vast field of study for the development of mathematical models in the area of systems biology¹⁷. The FKN model is extremely robust; however, it considers many chemical species as dynamic variables of the system, which makes it difficult to handle from an analytical point of view^15,16. To solve this problem, in 1974, Field and Noyes proposed a much simpler mechanism to describe reaction dynamics based on the FKN model, also known as Oregonator because it was created at the University of Oregon. However, said mathematical model was constructed from a reaction mechanism that consists of only five irreversible chemical reactions and that describes the BZ reaction qualitatively¹⁸. Years later, in 1990, Györgyi, Turányi, and Field proposed a detailed mechanism of the BZ reaction, that consisted of 80 chemical reactions and 27 chemical species and which is currently the most accepted reaction mechanism¹⁹. The BZ mechanism has two subsets, the first consists of a set of inorganic chemical reactions while the second one consists of organic chemical reactions. Just one year later, in his seminal work A history of chemical oscillations and waves, Zhabotinsky presented a mathematical model in differential equations based on a subset of chemical reactions from the work of Györgyi et al., and mentioned that it is possible to explain BZ reaction dynamics with said subset¹². At this point, the following questions arise: How reliable and valid is to reduce an extremely complex reaction mechanism to a subset of reactions? What criteria should be used to make this reduction? Could reducing the mechanism eliminate important system information? Is it possible to use mathematical tools to carry out this reduction properly without losing information? And if this is possible, can these mathematical techniques be used to study biochemical regulatory mechanisms to identify the most critical system variables?

To answer these questions, it is necessary to under-stand that DST is not the only mathematical tool used to study the complex mechanisms of biochemical regulation²⁰. There are tools such as models of agents and cellular automata that have been used to explain at least qualitatively the phenomena that emerge from different physiological processes²⁰. However, the tool that has attracted the most attention today is the Complex Networks Theory (CNT) based on graph theory^21,22. The CNT has been mainly developed by physicists, with the pioneering works of Barabasi et al.,²³. This theory has allowed us to understand the emergence of extremely complex behaviors in systems that involve a large number of variables because, through statistics, it enables us to study the structural properties of the networks to identify the variables or entities with greater relevance in the system of interest²⁴. Multiple measures of centrality can be used to determine the level of importance of a variable in a complex network, such as degree, clustering, hub, authority, etc.,²⁵. In Biomedical research, CNT is widely used to study many phenomena including dis-ease propagation, genetic regulation networks, protein-protein interaction networks, and identification of possible therapeutic targets in complex biochemical reaction mechanisms²⁵. In their excellent work, Costa et al. describe the CNT as a vital tool for Systems Biology²⁶.

It is well known that the emergence of chronic degenerative diseases (such as insulin resistance, diabetes mellitus II, cancer, cardiovascular diseases among others) result from the mismatch of a wide variety of physiological processes that are in turn regulated by an intricate network of biochemical reactions, which makes it difficult to study the interaction, coupling, and activation that may exist between different physiological processes^27,28,29,30. The CNT could facilitate the study of different biochemical processes related to each other taking these processes as complex networks where the participating chemical species can be considered nodes or vertices, and all possible physicochemical interactions between them as links or edges. It is possible to identify the chemical species with greater relevance by determining network centrality properties, and applying the standard chemical kinetics techniques (CK)³¹ gives rise to mathematical models in differential equations that facilitate the study of coupling and synchronization phenomena between different physiological processes and their possible relationship with various pathologies (see Figure 1).

FIGURE 1 General diagram for obtaining a system of differential equations from a biochemical system. It is possible to build complex networks from biochemical reaction mechanisms and, through the network's structural properties, identify the most relevant variables to use CK standard tools to build a dynamic system that is representative of this process. Finally, coupling the systems of equations obtained employing a bidirectional or unidirectional coupling to study possible synchronization effects and their relationship with various pathologies. The block diagrams of the biochemical processes are only schematic representations of Insulin Resistance and Diabetes Mellitus II. The colors of the blocks and lines are also a graphic representation of the chemical species and the possible interactions present in these processes. However, they do not describe the actual biochemical process. For the construction of the network diagrams, the Gephi software was used (https://gephi.org/)³⁴. 

It is necessary to identify when a group of physiological processes is coupled or synchronized, which requires tools that permit identifying such phenomena. Thanks to the development of machine learning, it is possible to create models capable of learning to recognize or identify a series of patterns with high precision through the acquisition of experience (data)³². Machine learning can be divided in supervised and unsupervised learning and can be further classified into different combinations of these^32,33. In general, supervised learning consists of providing a series of input data with their respective label to a model so that the model is generalized and subsequently allows the label to be predicted knowing only the input data or characteristics^32,33. On the other hand, unsupervised learning consists of providing only the input data or characteristics to the model without its label, to identify specific patterns of information^32,33. Each supervised learning model must go through a training, validation, and testing process to ensure generalization^32,33. Supervised learning models can be divided into two large groups: classification models and regression models; in the former, the variable is to be predicted as a qualitative or categorical variable. In regression models, the variable to be predicted is a quantitative variable, either continuous or discrete^32,33. Hereby, the following question arises: Is it possible to use a supervised learning model to predict the coupling or synchronization of biochemical processes represented by systems of differential equations?

To answer this question, it is necessary to ask another one: How can data be provided to the supervised learning model to achieve the prediction of a coupling state? The answer is that by modeling biochemical processes as systems of differential equations and coupling them unidirectionally or bidirectionally, it is possible to numerically solve these models, from which the change in the concentration of chemical species with respect to time is obtained. Therefore, it is possible to determine if, for any value of the coupling parameter, these chemical species are synchronized. The degree of synchronization can be determined by evaluating some nonlinear metric, such as the fractal dimension of the time series, evaluation of Lyapunov exponents, entropy, the study of the synchronization variety in the phase space of the variables under study, etc.^35,36,37,38, then it is possible to use some of these metrics to determine, through a supervised learning model, whether these systems are coupled or not. It is also possible to obtain images from the numerical solution of the systems of differential equations; such is the case of recurrence diagrams^39,40 or Gramian angular field images (gaf)⁴⁰, which can be used with a supervised learning model to deter-mine if a group of variables, in this case, chemical species related to important biochemical processes, are coupled or not, and thereafter, to associate a possible state of synchronization with various pathologies. Currently, it is possible to monitor blood concentration of a variety of chemical species related to critical metabolic processes, which are associated with different pathologies^41,42,43,44, making it possible to construct time series and use them to build a database of gaf images to train, validate and test a supervised learning model to identify possible synchronization states and their relationship with pathologies; however, clinical trials can be expensive and, as a result, conducting this type of study is not practical^41,42,43,44. On the other hand, using mathematical models that represent the dynamics of these variables and using their numerical solution to obtain their respective gaf images to train, validate and test the supervised learning model, offers a solution to this problem. The ability of each mathematical model to represent a biochemical process depends on its degree of complexity, i.e., whether the model involves not only metabolic biochemical processes but also epigenetic processes⁴⁵.

To analyze the possibility of using CNT in biochemical regulatory mechanisms for identifying the most critical variables of the system and building a system of differential equations that model it, the BZ reaction was used as a study model since it is currently the most complex oscillating reaction discovered so far and it has been used to study the biochemical feedback mechanisms present in the Krebs cycle^{12,13,14,15,16,17,18,19}. On the other hand, to emulate the effects of coupling and synchronization, the obtained model was bidirectionally coupled with other models associated with oscillatory biochemical processes, and its numerical solution was used to obtain the gaf images. Finally, a supervised learning tool was used to identify the type of models coupled using the gaf images.

Materials and methods

For this study, the general mechanism of the BZ reaction proposed by Györgyi et al.¹⁹ was used. A network was built considering two types of nodes or vertices, the "chemical species" type nodes, and the "chemical reaction" type nodes. The link or edge between species and chemical reactions is given by the reaction rate constant of each of the reactions, which gives place to a "bipartite" network^{20,21,22,23,24,25}. Once the net-work was built, the structural property authority was evaluated using the algorithm proposed by Kleinberg, which was initially proposed to determine the level of importance and information flow of websites, to reveal the sites with the highest traffic in a virtual hyperlink environment⁴⁶. In this work, the structural property authority was used to determine the importance of each of the chemical species involved in the BZ reaction mechanism. To build a system of nonlinear differential equations capable of describing the BZ reaction mechanism, it was assumed that the most relevant variables have the most significant flow of "chemical" information. Once the most critical chemical species in the reaction mechanism were determined under the authority criteria, the mathematical model was constructed using standard techniques of CK³¹. When the model was obtained, the effects of synchronization that could emerge due to its coupling with different chemical oscillators were studied to emulate the synchronization processes present in different biochemical systems in mammals⁴⁷. To this end, a bidirectional coupling was carried out with three different models⁴⁸. The first one was an identical model, the second one was the model proposed by Levefer also called Brusselator^49,50, which shows interesting autocatalytic processes and, lastly, the model proposed by Selkov that describes the oscillatory behaviors present in glycolysis^51,52. For each case, the coupling parameter was varied and the numerical solution of the systems of coupled differential equations was obtained using the standard fourth-order Runge Kutta technique⁵³. Once the time series of the numerical solution of the differential equations were obtained, the image of the time series of the bidirectionally coupled variables was constructed using the gaf technique⁵⁴, which was described extensively by Wang et al⁵⁵. In general, this technique consists in representing a time series in a polar coordinate system, each X_i value of the time series X is scaled so that all values are in the interval [-1,1] or [0,1], subsequently, with the rescaled values the angular cosine is obtained φi=arccos⁡xi,-1≤xi≤1,xi∈X~ and the coordinate rr=tiN,ti∈N, where t_i corresponds to the time interval between each value of the time series and N is a constant factor⁵⁵. Lastly, the addition/subtraction between each point is determined to identify possible temporal correlations within different time intervals⁵⁵. It is possible to obtain two types of gaf, the sum (gasf= cos (φ_i +φ_j)) and the subtraction (gadf= sin (φ_i +φ_j)). These values are used to construct the Grammar matrix that is used to obtain the image (see Figure 2)⁵⁵. This technique in combination with supervised learning tools has been used to study time series of electroencephalograms⁵⁶, electrocardiograms⁵⁷, signals obtained from biosensors⁵⁸, etc.

FIGURE 2 General diagram for obtaining the gaf images. The time series are rescaled, later these are represented in a polar coordinate system and finally, the gaf images are obtained. 

After the gasf and gadf images were obtained, the supervised learning model was generated using the Orange Data Mining software (https://orange.biolab.si/)^58,59, which is a visual toolbox where it is possible to build workflows that allow the use of widgets for analysis and data processing, supervised and unsupervised learning tools, data visualization and model evaluation. It also has extensions for text mining, spectroscopy, complex network analysis, time series, bioinformatics, and image analysis^59,60,61. The transfer learning technique, which is an artificial intelligence technique that consists of pre-training a model with an extensive database and the experience gained from said training to apply it to another problem that may be completely different, was used to process the gaf images^61,62,63. This technique is used in image processing as follows: deep convolutional neural networks (CNN) are used, which are pre-trained with a large number of images of all kinds, later, activations of the penultimate layer of the model (CNN codes) are used to represent the images with vectors (embedded), i.e., CNN is used as a feature extractor or descriptor, allowing supervised learning models to be used and thus obtaining high precision values in image classification^61,62,63. In Orange, it is possible to embed images using different CNNs, including Google's Inception V3 neural network that has been pre-trained with the ImageNet database consisting of 1.2 million images. The neural network has 2048 nodes in its penultimate layer, so each image rep-resents it with a vector of dimension 2048^33,64 (CNN by default in Orange). On the other hand, embedding images can also be done with CNN SqueezeNet, which is a much simpler network than Inception V3; nevertheless, it achieves a precision close to that of CNN AlexNet on the ImageNet database, this CNN represents each image as a vector of dimension 1000⁶⁵. It is also possible to embed with CNN's VGG16 and VGG19 proposed by the Visual Geometry Laboratory of the University of Oxford⁶⁶, in the same way, pre-trained with the ImageNet database, the CNN DeepLoc pre-trained with 21882 images⁶⁷ and CNN Painters, a pre-trained network with 79433 images^59,60,61. For embedding the images, Orange sends them to an external server, except for CNN SquezeeNet, which is done locally, making embedding faster. In this work, the gaf images were embedded using CNN Google Inception V3 and SquezeeNet. Subsequently, a supervised learning model was trained as a classifier of the type of coupled oscillators using the gaf images. A logistic regression with Ridge penalty or regularization (L2 = 1) was used as a classification method. Six different evaluation techniques were used to train, validate and test the classifier: 1) the standard stratified holdout technique (70% training set / 30% test set), the training and testing subsets are repeated 10 times (A), 2) a 3-fold stratified cross-validation ( B), 3) a 5-fold stratified cross-validation (C), 4) a 10-fold stratified cross-validation (D), 5) a 20-fold stratified cross-validation (E) and 6) the leave-one-out cross-validation (F), to identify the best evaluation technique and avoid overfitting^33,68,69. The previous procedure was repeated implementing a principal component analysis (PCA) using a number of components such that they explain 95% of the total variance of the images (21 principal components for Google Inception V3 and nine principal components for SquezeeNet). After embedding the images with the CNN, the same six evaluation techniques were used to study the effect of reducing dimensions and to avoid overfitting^33,68,69. For both procedures the classifier confusion matrix was obtained to determine the classification metrics: Classification Accuracy (CA), Precision (P), Recall (R), and F1-score (F1)^33,70.

Results and discussion

Figure 3a is a diagram that represents the complex network for the general mechanism of the BZ reaction proposed by Györgi et al. The nodes represent each of the 27 chemical species and 80 chemical reactions participating in the reaction mechanism, and the links represent the reaction rate constant magnitude. The largest nodes, which are shown in Table 1, represent the nodes with the highest numerical value of authority. Figure 3b shows the network authority property histogram. Of the 107 nodes, 49 (i.e. 45.79%) have an authority value between 0 and 0.1 (none being actually zero), 40 have an authority value between 0.1 and 0.2 (37.38%), 9 nodes between 0.2 and 0.3 (8.41%), 3 nodes between 0.3-0.4 (2.8%), 3 nodes between 0.4 and 0.5 (2.8%), 2 nodes between 0.5 and 0.6 (1.87%), and one between 0.9 and 1 (0.93%). This last one is the most important node, considering the authority centrality criterion. In the histogram, it is easy to distinguish that there are few nodes with a high authority value, making it possible to identify the chemical species with greater importance, possibly due to power-law-like behavior⁷¹. Power laws are present in a wide variety of biological phenomena and are related to processes of a universal nature⁷².

FIGURE 3 a) The complex network of the chemical mechanism of Györgyi et al¹⁹. The larger nodes correspond to those that have greater authority value; b) Histogram of structural property authority. It can be observed that most of the nodes have an authority value below 0.3 (91.58%), however, there are few nodes with authority values above 0.3 (8.42%), which can be considered as the most important nodes under the criterion of centrality authority. Gephi software was used for the construction of the network diagram (https://gephi.org/)³⁴. 

As expected, the chemical node or species with higher authority is the H+ cation because the chemical reaction at the experimental level is carried out in the presence of sulfuric acid (H₂SO₄), which is a source of H+ and allows the formation of bromous acid (HBrO₂), an essential variable for the chemical feedback mechanism that enables the existence of periodic behaviors^12,13,14. Furthermore, we can note that H₂O is the second chemical species with higher authority since the precursor solutions for the chemical reaction use deionized water as a solvent¹². The third most crucial chemical species is the carboxy radical (*COOH) inasmuch as the breakdown of malonic acid as a precursor agent results in short-lived chemical species^12,19. The presence of Br* radicals is mainly because during the progress of the reaction molecular bromine (Br₂) occurs due to the bromate anion precursor agent (BrO₃), said Br₂ results in the formation of bromide ions (Br-) and subsequently to Br* radicals, which are fundamental chemical species for the feedback mechanism^12,19. Because the chemical reaction includes organic components, specifically malonic acid and its derivatives, it has carboxyl or carbonic acid groups that can be decomposed into carbon dioxide (CO₂), which is seen in the form of gas bubbles during the reaction^12,13,14.

Lastly, we can see that the oxidized form of the Cerium metal (Ce⁴⁺) is also among the ten nodes with the highest authority value, this results from said chemical species acting as a catalytic agent in the chemical reaction and that alongside its reduced form (Ce³⁺), it gives rise to the typical color changes of the BZ reaction^12,18. With these results, it is possible to observe that applying CNT to the BZ reaction mechanism proposed by Györgyi et al.¹⁹ permits identifying the chemical species with greater relevance under the criterion of authority. Such chemical species appear naturally in the mechanism described by Zhabotinsky, either in its ionic form or in the form of a precursor compound¹², so it is possible to approximate the BZ reaction to said subset of reactions. Still, it is essential to mention that it is possible to use other centrality criteria to identify the nodes with greater relevance²⁵.

The mechanism proposed by Zhabotisnky uses Iron (Fe) as a catalyst instead of Cerium (Ce)¹². Below is the chemical mechanism studied by Zhabotinsky replacing Fe with Ce.

Equations 8

εdu1dt=u11-u1+fw1+αq-u1q+u1+ku2-u1,

dw1dt=u1-w1.

du2dt=v-u2w2γw2γu2+w2γ+1+ku1-u2,

dw2dt=βu2w2γw2γu2+w2γ+1-ηw2.

In Figure 4a, the numerical solution of the system of Equations 6 is shown. A full synchronization can be seen between variables u₁ and u₂, i.e., the oscillation rhythms are completely coupled^48,75. This type of behavior is common in biochemical processes carried out by cells that share the same microenvironment and can be stimulated either by a chemical or physical mechanism^76,77. On the other hand, in Figure 4b the gasf image of the variable u₁ can be seen, while in Figure 4c, the gadf image of the variable u₂ is shown. A characteristic pattern of bidirectional coupling between identical oscillators with periodic behavior can be seen in both figures.

FIGURE 4 a) Numerical solution of the system of Equations 6 with ε₁= ε₂= 0.3, α₁= α₂= 0.03, q₁= q₂= 0.015, ƒ₁= ƒ₂= 1.0 and k= 1.0. With u₁ (0)= 0.3, w₁ (0)= 0.5, u₂ (0)= 0.8 and w₂ (0)= 0.7, tƒ=1 50 units of dimensionless time and dt=0.001. An oscillatory behavior is shown for the 4 dynamic variables of the system and full synchronization is observed; b) gasf image of the variable u₁; c) gadf image of the variable u₂. For both images, a characteristic pattern of periodic systems can be seen. 

In Figure 5a, the numerical solution for the system of Equations 7 is shown. Almost complete synchronization can be seen between the variables u₁, w₁ and u₂^48,75. Anyway, this system of coupled differential equations is very sensitive to the value of the coupling parameter. Also, Figure 5b shows the gasf image of the variable u₁, while Figure 5c shows the gadf image of the variable u₂. In both figures, a characteristic pattern for a quasi-periodic system can be seen⁷⁸.

FIGURE 5 a) Numerical solution of the system of Equations 7 with ε= 0.3, α= 0.03, q= 0.015, ƒ= 1.0, a= 1, b= 2.5 and k= 2.5. With u₁ (0)= 0.3, w₁ (0)= 0.5, u₂ (0)= 0.8 and w₂ (0)= 0.7, tƒ= 150 units of dimensionless time and dƒ=0.001. A quasi-periodic behavior is shown for the 4 dynamic variables of the system and almost complete synchronization is observed between the variables u₁, u₂ and w₁; b) gasf image of the variable u₁; c) gadf image of the variable u₂. The images reflect the quasi-periodic behavior of the variables u₁ and u₂. 

In Figure 6a, the numerical solution of the system of Equations 8 is shown. The phase synchronization between the variables u₁, u₂ and w₁^48,75 can be seen. The NOM model oscillation rhythm causes variable u₂ to enter in-phase synchronization; however, the variable w₂ decays completely to zero. In the Selkov model, the variable u₂ is related to the concentration of ATP (adenosine triphosphate), while the variable w₂ is related to the concentration of ADP (adenosine diphosphate)^51,52, so when said system is coupled with the NOM model, the synchronization process that could exist between the biochemical mechanism of glycolysis and the Krebs cycle is emulated.

FIGURE 6  Figure 6. a) Numerical solution of the system of Equations 8 with ε= 0.3, α= 0.03, q= 0.015, f= 1.0, ν= 0.0285, η= 0.1, β= 1.0, γ= 2 and k= 0.1. With u₁ (0)= 0.3, w₁(0)= 0.5, u₂ (0)= 0.8 and w₂ (0)= 0.7, tƒ= 150 units of dimensionless time and dƒ= 0.001. Oscillatory behavior is distinguished for the variables u₁, u₂, and w₁, as well as phase synchronization between them. It is recognized that the variable w₂, for this case of the value of the parameters, tends to zero; b) gasf image of the variable u₁; c) gadf image of the variable u₂. The gasf image shows a typical pattern of a periodic oscillatory system. On the other hand, although the variable u₂ shows an oscillatory behavior, it is also coupled with the variable w₂, which is why a well-defined pattern in the gadf image is not clearly distinguished. 

Moreover, in Figure 6b the gasf image of the variable u₁ can be seen, while in Figure 6c the gadf image of the variable u₂ is shown. It can be distinguished that the gasf field shows the typical oscillatory behavior, while the image of the gadf field does not clearly show a distinct pattern.

The numerical solution of the equations without coupling can be found in the supplementary material, as well as their respective gasf and gadf images (see Figures S1, S2 and S3). The coupling of chemical oscillators to emulate the synchronization in biochemical processes, is of vital importance since it allows us to understand the complex dynamics that underlie the feedback processes present in biological systems. Mismatches in the oscillation rhythms of physiological processes can lead to a wide variety of metabolic problems⁷, so it is essential to know the type of coupling that can exist between oscillators (unidirectional, bidirectional, linear, nonlinear, etc.,⁴⁸.) and whether these are coupled or not at any given time.

To determine coupled oscillators type, the coupling parameter was varied between one and ten for the system of Equations 6, between one and ten for the system of Equations 7 and lastly, between 0.1 and 1 for the system of Equations 8 (because this system is highly sensitive to bidirectional coupling), to generate the images of the time series using the gaf technique. Once obtained, the images were embedded to subsequently train the supervised learning model using the procedure described in the materials and methods section (see Figure 7). Figure 8 shows the general Orange workflow for this procedure. In the supplementary material (Tables TS-1 and TS-2), the results of the evaluation of different models of supervised learning are shown (Multilayer Perceptron Neural Network (MLP) (Hyperparameters: Optimizer: L-BFGS-B, activation function: Logistics, regularization L2= 1, neurons in one hidden layer: 10, tol= 1E-4, max_iter= 200), K-Nearest-Neighbors (Hyperparameters: Number of nearby neighbors: 5, metric: Euclidean, wight= uni-form) and Support Vector Machine (SVM) (Hyperparameters: Cost= 1.0, kernel: RBF, tol= 1E-3, iteration limit= 100))^33,68, which were used to classify the type of coupled oscillators using the gaf images, to compare them with those obtained with the logistic regression model. When training the different super-vised learning models using the vectors provided by CNN SquezeeNet as image descriptors, the logistic regression model presented the highest CA and F1 val-ues for all evaluation techniques, obtaining the maxi-mum values in the methods A, D, E and F (CA= 0.983 and F1= 0.983). In contrast, the minimum values were obtained from methods B and C (CA= 0.966 and F1= 0.967) (see supplementary Tables ST-1). When using CNN Inception V3 to embed the images, similar values of CA and F1 were obtained for the logistic regression models, MLP and KNN for method A (CA= 0.961, F1= 0.961), methods B, C, D, E and F (CA= 0.950 and F1= 0.949-0.950). The SVM model presents the highest CA and F1 for methods B and C (CA= 0.966 and F1= 0.966). Notwithstanding, for methods D, E and F, the SVM model presents similar CA and F1 than the rest of the models (CA= 0.950 and F1= 0.949) (see supplementary Tables ST-1). Likewise, it is possible to observe that there are no major differences between the results obtained by CNN's Inception V3 and SquezeeNet. Godec et al., use transfer learning to embed biomedical images and mention that this technique allows them to obtain high precision values in classification models using small databases⁶¹. Godec et al., also found no major differences in the CA and F1 values of the logistic regression classifier when using either the CNN Inception V3 or the CNN SquezeeNet to embed the images⁶¹.

FIGURE 7 General procedure to obtain the classification model of the type of coupled oscillators using the gaf images. a)The numerical solution of the coupled oscillator systems is obtained by setting the parameters of the models and varying the coupling parameter, b) subsequently gaf images are obtained from the time series of the numerical solution of the coupled oscillators, c) after building the database, then embedding the images using Google's CNN Inception V3 or SquezeeNet and using PCA for dimension reduction (optional). Lastly, d) the classification model is evaluated using the techniques: stratified holdout (70% training set / 30% testing set), 3, 5, 10, 20-fold stratified cross-validation, and leave-one-out cross-validation. 

FIGURE 8 Orange workflow example for training the logistic regression model as a gaf image classifier. a)In the Orange interface, the images are imported first from a local folder using the "Import images" widget, then the images are embedded using the "Image embedding" widget and finally the training, validation and supervised learning model test. b) Window that shows the "Image Embedding" widget, which shows the CNNs that can be detected to extract the characteristics or descriptors of the images. The CNN of Google Inception V3 is shown by default. c) Window that shows the "Logistic Regression" widget. It is appreciated that it is possible to use the regularization of the sea Lasso (L1) or Ridge (L2). d) Window that shows the widget "Test and score." This widget shows the different evaluation techniques, which are: "k-fold cross-validation", "k-fold stratified cross-validation", "Random sampling", and "Leave-one-out cross-validation". Also, it shows the classification metric by class or average. e) Window that shows the "Confusion matrix" widget, which allows observing the confusion matrix of the supervised learning model. The widget allows displaying the confusion matrix in percentages or the number of images classified correctly and incorrectly. 

It is worth mentioning that sometimes when the data-base used for training supervised learning models has a higher number of characteristics or descriptors compared to the sample size (as is often the case in biomedical databases), or is unbalanced, i.e., there is a greater amount of data from one class than from the rest, it is possible to overfit the model, leading to erroneous results^33,61,68,69. However, it is possible to prevent overfitting essentially through two procedures, the first is to use more data for training, decrease statistical bias and decrease the number of characteristics or descriptors, the second one is to limit the complexity of the model of supervised learning, employ regularization, either penalty L2 (Ridge), L1 (Lasso) or ElasticNet (L1 and L2 simultaneously) and use assessment techniques such as stratified cross-validation or leave-one-out cross-validation^{33,68,79,80,81,82,83}. Therefore, we have also implemented a PCA after embedding the images with CNN's to study the effect of the reduction of dimensions in the classification of the type of coupled oscillators, using the same models of supervised learning and the same evaluation methods.

When performing the PCA implementation after embedding the images with CNN SquezeeNet using nine main components, which explain 95% of the total variance, the logistic regression model obtains the highest CA and F1 value for method A (CA= 0.983 and F1= 0.983) nevertheless, it also shows the best classification metrics for methods C, D, E and F (CA= 0.966 and F1= 0.967), while for method B, the SVM model presents the best classification metrics with a CA= 0.966 and F1= 0.966 (see supplementary material Tables ST-2). Moreover, when a PCA is implemented after embedding with CNN Inception V3 using 21 main components, which explain 95% of the total variance, the logistic regression model shows the best classification metrics for the evaluation methods A, C, D, E and F, obtaining a value of CA= 0.961 and F1= 0.960 for method A, while for methods C, D, E and F, a CA= 0.950 and F1= 0.949 (see tables of supplementary material ST-2) were obtained. Regardless, the SVM model presents the highest values of CA and F1 for method B (CA= 0.966 and F1= 0.966). Having said that, it is possible to note that there are no major differences in the classification metrics when using either the CNN SquezeeNet or the CNN Inception V3 to embed the images and, at the same time, there are no major differences in the metrics of classification when using the descriptors connected directly from CNN or when the main components are extracted from them. However, the computational time required for the evaluation of the models when using PCA is shorter. The supervised learning models present similar values in the classification metrics; nevertheless, the logistic regression model has the least complexity because it only uses one hyperparameter, which is used as a regularization or penalization^33,67,84,85. Therefore, logistic regression can be used as a classification model for the type of coupled oscillators. Table 2 shows the classification metrics for the logistic regression model for each of the evaluation methods, using descriptors obtained directly from CNN's as training data.

TABLE 2 Classification metrics for the logistic regression model using the characteristics obtained directly from the CNN for the evaluation. 

Evaluation Technique	CNN
	Inception V3				SquezeeNet
A	CA	F1	P	R	CA	F1	P	R
B	0.961	0.961	0.963	0.961	0.983	0.983	0.984	0.983
C	0.950	0.949	0.953	0.950	0.966	0.967	0.969	0.966
D	0.950	0.949	0.953	0.950	0.966	0.967	0.969	0.966
E	0.950	0.949	0.953	0.950	0.983	0.983	0.984	0.983
F	0.950	0.949	0.953	0.950	0.983	0.983	0.984	0.983
A	0.950	0.949	0.953	0.950	0.983	0.983	0.984	0.983

Likewise, Table 3 shows the classification metrics for the same model, using the descriptors obtained from the PCA as training data. At this point, it is natural to ask what evaluation method should be used if all methods have similar ranking metrics. For this work, we chose method D because multiple experiments have been carried out that demonstrate that the best way to obtain high values in the metrics, be it classification or regression, is using stratified 10-fold cross-validation, even when there is the possibility of computation to increase fold number in the evaluation of supervised learning models^33,86. In addition, as can be seen in Tables 2 and 3, for method D the same values of the classification metrics are obtained using the descriptors extracted directly from CNN Inception V3 and those obtained from the implemen-tation of the PCA ( CA= 950, F1= 0.949-0.950) for train-ing. When using the descriptors extracted directly from CNN SquezeeNet for training, the classification metrics CA= 0.983, and F1= 983 were obtained, while those obtained due to the implementation of the PCA are CA= 0.967 and F1= 0.967, which means that there is no significant difference. In conclusion, there are no major differences between the use of the descriptors extracted from the implementation of the PCA after embedding the images with one or the other CNN. This shows that the dimensions reduction does not substantially affect the precision of the supervised learning model and, conversely, allows for a better generalization of it⁶¹.

TABLE 3 Classification metrics of the logistic regression model using the characteristics obtained from the application of the PCA. 

Evaluation Technique	CNN
	Inception V3				SquezeeNet
	(21 PC, explained variance: 95%)				(9 PC, explained variance: 95%)
A	CA	F1	P	R	CA	F1	P	R
B	0.961	0.961	0.963	0.961	0.983	0.983	0.984	0.983
C	0.950	0.950	0.953	0.950	0.950	0.951	0.957	0.950
D	0.950	0.950	0.953	0.950	0.967	0.967	0.970	0.967
E	0.950	0.950	0.953	0.950	0.967	0.967	0.970	0.967
F	0.950	0.950	0.953	0.950	0.967	0.967	0.970	0.967
A	0.950	0.950	0.953	0.950	0.967	0.967	0.970	0.967

Figure 9 shows the confusion matrix of the logistic regression model trained with the descriptors obtained from applying the PCA after being embedded with CNN Inception V3 and using the evaluation method D. 18 of the 20 images of the coupling of the NOM and Brusselator Oscillators (O-B) have been correctly classified, while two images have been erroneously classified, one as a coupling between the NOM and Silkov (O-G) oscillators and another as a coupling between the identical NOM (O-O) oscillators. Furthermore, 19 of the O-G coupling images have been correctly classified, while one has been erroneously classified as O-O coupling. All O-O coupling images have been correctly classified.

FIGURE 9 Confusion matrix of the logistic regression method using the CNN Inception V3 as an image descriptor extractor and applying a reduction of dimensions using the PCA technique (21 PC, explained variance: 95%) (evaluation method D). It is observed that 18 of the 20 images of the O-B coupling have been classified correctly, while two have been erroneously classified, one as the O-G coupling and the other as the O-O coupling. On the other hand, the coupling shows that 19 of the 20 images of the O-G coupling have been correctly classified, while one has been classified as an O-O coupling. All images in the O-O coupling have been correctly classified. 

Figure 10 shows the confusion matrix of the same model trained with the descriptors obtained from applying the PCA after being embedded with CNN SquezeeNet and using the evaluation method D. In the confusion matrix, it is possible to observe that 19 of the 20 images of the O-B coupling have been correctly classified, while one has been incorrectly classified as O-O coupling. Similarly, 1 of the 20 images of the O-G coupling has been incorrectly classified as O-O coupling. All images in the O-O coupling have been correctly classified. The decision to use one or the other CNN for embedding the images will depend on whether, as users, we want our images to be sent to an external server for embedding. For privacy and security rea-sons, we prefer them to be embedded locally⁶¹.

FIGURE 10 Confusion matrix of the logistic regression method using as descriptor extractor of CNN SquezeeNet images and applying a reduction of dimensions using the PCA technique (9 PC, explained variance: 95 %). It is observed that 19 of the 20 images of the O-B coupling have been correctly classified, while one has been classified as O-O coupling. While, 1 of the 20 images of the O-G coupling has been incorrectly classified as O-O coupling. All images in the O-O coupling have been classified correctly. 

Conclusions

In the framework of Systems Biology, mathematical modeling of biochemical mechanisms involved in different physiological processes is of vital importance because it allows us to understand the non-linear dynamics that underlie these phenomena. This is why the use of mathematical tools and computational systems for the analysis of the complex feedback mechanisms present in living systems is necessary. The CNT is a mathematical tool that allows studying these mechanisms with a holistic approach and provides valuable information on each of the entities that make up the system^26,87.

When determining the authority structural property of the complex network obtained from the BZ reaction mechanism proposed by Györgyi et al., and using it as the centrality criterion, the variables with the highest relevance were identified, i.e., those chemical species that have the greatest flow of information and that could participate in the emergence of collective properties of the system. Identification of these variables led to the construction of a nonlinear system of differential equations similar to the reduction of the FKN model proposed by Field and Noyes (Oregonator) and which also explains the phenomenology of the BZ reaction. Hence, this result answers the question of using mathematical tools to reduce complex reaction mechanisms without losing generality. Therefore, it is possible to use this methodology in the study of nonlinear dynamics present in biochemical and physiological processes.

On top of that, by applying this methodology to biological systems, it is possible to translate any biochemical or physiological process to a mathematical model and study the phenomena of synchronization between different regulatory mechanisms⁸⁸ to decipher the complex dynamics that underlie living systems with a systemic approach.

The effect of coupling between oscillators of different nature can be clearly seen in the images obtained using the gaf technique, which can be used to train a supervised learning model to classify the type of coupled oscillators. The extraction of descriptors from gaf images through pre-trained CNNs (transfer learning) allows obtaining high precision values in the evaluation of different classification models; however, it is also possible to couple the pre-trained CNNs with the PCA to obtain high values in the classification metrics, comparable with the values of the metrics obtained by using only pre-trained CNNs as a descriptor extraction method. In particular, using the CNN's Inception V3 and SquezeeNet as extractors for descriptors of gaf images and obtaining the principal components of these descriptors, allows training classification models such as logistic regression and obtaining CA and F1-score values above 0.94 for different evaluation methods.

All things considered, the methodology proposed in this work can facilitate the determination of synchronization and desynchronization states in complex real biochemical and physiological mechanisms to recognize a possible correlation between these states and the emergence of different complex diseases.