GENAVOS: A New Tool for Modelling and Analyzing Cancer Gene Regulatory Networks Using Delayed Nonlinear Variable Order Fractional System

: Gene regulatory networks (GRN) are one of the etiologies associated with cancer. Their dysregulation can be associated with cancer formation and asymmetric cellular functions in cancer stem cells, leading to disease persistence and resistance to treatment. Systems that model the complex dynamics of these networks along with adapting to partially known real omics data are closer to reality and may be useful to understand the mechanisms underlying neoplastic phenomena. In this paper, for the ﬁrst time, modelling of GRNs is performed using delayed nonlinear variable order fractional (VOF) systems in the state space by a new tool called GENAVOS. Although the tool uses gene expression time series data to identify and optimize system parameters, it also models possible epigenetic signals, and the results show that the nonlinear VOF systems have very good ﬂexibility in adapting to real data. We found that GRNs in cancer cells actually have a larger delay parameter than in normal cells. It is also possible to create weak chaotic, periodic, and quasi-periodic oscillations by changing the parameters. Chaos can be associated with the onset of cancer. Our ﬁndings indicate a profound effect of time-varying orders on these networks, which may be related to a type of cellular epigenetic memory. By changing the delay parameter and the variable order functions (possible epigenetics signals) for a normal cell system, its behaviour becomes quite similar to the behaviour of a cancer cell. This work conﬁrms the effective role of the miR-17-92 cluster as an epigenetic factor in the cancer cell cycle.


Introduction
The gene regulatory network (GRN) is a set of molecular regulators that interact with each other and with other cellular materials to control the expression level of mRNA and proteins. The regulator can be DNA, RNA, proteins, and their complexes. This interference can be direct or indirect (via transcribed RNA or translated protein) [1]. In the case of genomics, the concept of genes has been developed into two categories: protein-encoding genes and non-encoding genes. Micro-RNAs (miRNA) and long non-coding RNAs (LNC-RNA) are non-coding genes. By this latter definition, genes are regulated by other genes, proteins, metabolites, and DNA, as well as epigenetic factors in a gene regulatory network [2]. Epigenetics are defined as mechanisms that lead to heritable changes in gene functions without alteration in the DNA sequence [3]. Among these mechanisms involved in GRNs are DNA (cytosine) methylation and hydroxymethylation, chromatin remodelling, and non-coding RNAs [4]. Also, environmental factors, such as drugs and chemicals, temperature, light, and circadian rhythm, affect GRNs [5]. Essentially, cancer is abnormal process of alteration in the cell molecular regulation network [6]. The role of GRNs can go beyond cancer pathogenesis to cell invasion, metastasis, and drug resistance. These key roles are consistent with the cancer stem cells (CSC) theory [7]. CSC theory states that a small number of cancer cells act as stem cells that reproduce and maintain cancer [7]. This theory suggests that tumorigenesis may be regulated by maintaining a balance between asymmetric and symmetrical cell division. Thus, mutations affecting this balance can lead to abnormal expansion of stem cells [8]. Disruption in GRNs regulates asymmetric cell division in CSCs through genetic factors (such as genes and proteins) and epigenetic factors (such as specific miRNAs and long non-coding RNAs) [7,8]. When asymmetric division occurs in many cell divisions, a population of cancer cells shows tumour heterogeneity, where many cells have slightly different genetic and phenotypic states. This is one of the key features that makes deep cancer adaptation possible [9]. Therefore, modelling GRNs is one of the important goals of system biology to study the creation, development, and treatment of cancer. Although an interaction between encoding genes occurs through proteins, or their metabolites, or compounds as well as epigenetic factors, reducing the GRN topology to a network of interactions between genes without considering intermediates between them is a simplified but logical approach to model their time-dependent dynamics [10,11]. An important reason for such simplification is the lack of information other than high throughput transcriptome level data, such as gene expression data obtained with the popular microarray or next generation sequencing (NGS) technology [10]. For this reason, gene expression time series data are commonly used to model and infer GRNs dynamics as well as dynamic biological processes [11][12][13]. Other omics data, such as epigenomics and proteomics, are commonly used to infer and predict the structure and topology of integrated GRN networks, including genetic factors (such as genes and proteins) and epigenetics, rather than modelling the dynamics of the system [3,13]. Although gene expression time series data are ultimately valid for describing phenomena observed with the transcription profile, they contain hidden and/or visible information from time-dependent dynamic biological processes such as the cancer cell cycle as well as latent information from other biological agents such as epigenetic factors, tumour interactions with the immune system, cell migration, etc. [14]. In other words, the gene expression signal in a gene regulatory network either affects and/or is affected by these factors [14,15]. Therefore, it is logical that if the gene expression time series data are used to model GRNs, the network topology should only include gene-gene interactions, but at the same time the model should justify phenomena beyond the transcriptional level as much as possible.
The first step in GRNs modelling is to determine the structure of the network. Then, a complete analysis of the dynamics of its components requires solving the Master Equations. By writing the master equations, this network is described as a dynamic system with nonlinear differential equations in state space. The state variable in GRN system can be time-varying concentration of mRNA (or gene expression) in a biological process such as a cell cycle [12]. In systems theory, GRNs are complex nonlinear dynamic systems with time delays, and complex nonlinear dynamics can be found even within the simplified and constructive infrastructure of these systems [16].
Delays are inherent in GRN systems and can lead to complex behaviours [16,17]. Transcription and translation are processes that take a certain amount of time to complete. Therefore, these two gene expression processes cause time delays in biochemical systems [18]. The nonlinearity of biochemical phenomena such as phosphorylation and dimerization, auto-regulation or feedback loops are known in GRNs [19]. Also, experimental studies have shown that the function of transcription has a sigmoidal form [20,21]. Although the listed intrinsic biological properties of GRN systems can be described by models of integer order deterministic or stochastic differential equations, GRN systems have other intrinsic biological properties that these models cannot describe. We explain it below.
It has been shown that transcription dynamics are slow and epigenetic dynamics such as DNA methylation dynamics are even slower [22,23]. Epigenetic mechanism by chromatin modification and/or DNA methylation may provide cellular biochemical memory by blocking or allowing transcription [24,25]. An essential role of epigenetic memory is that the state of gene expression must be maintained or re-established through each cell cycle. Likewise, an epigenetic state must be maintained through the phases of the cell cycle [26]. Another epigenetic control mechanism is through non-coding RNAs. The complete GRN model is computationally difficult to simulate and interpret its dynamic behaviour conceptually. Even if such a complete model exists, it will be difficult to validate it with real data because, as we said before, there is no ChIP-seq (chromatin immunoprecipitation) or ChIP-chip time series data as well as non-coding RNA expression time series similar to what is available for gene expression [23].
Fractional order systems are very suitable for modelling slow temporal dynamics of transcription and epigenetic. Also, fractional differentials can reflect the memory properties of GRNs. The fractional order system also provides a large degree of freedom, thus dealing with the nonlinear effects of GRNs. On the other hand, the dynamics of GRNs depend on the fractional order of state space equations. Although some optimal dynamic behaviours are obtained by selecting the appropriate fractional order, due to the high complexity of GRNs, there is a need for more flexible dynamics than fractional order equations [27]. This highly flexible dynamic can be obtained by generalizing fractional order systems to variable order fractional systems (VOFS). In the VOF-GRN system, the fractional order of the differential equations can depend on time, environmental and epigenetic factors. Variable order fractional derivatives/integrals can be operators to describe a system with time-varying memory [28]. These systems are more compatible with biological phenomena of epigenetics [29]. All of the above suggests that VOF systems provide a better description of the dynamics of GRNs. The main idea of this paper in using the VOF system to model GRNs is this hypothesis: Epigenetic signals that control the regulation of gene expression and create memory in the GRN system can be the variable order of the system as a separate input signal from Gene-gene networks but affecting them. It can also be seen that our method has an approach to approximate the time series of non-coding RNAs expression. Figure 1 shows the simplified biological system of the proposed method by simplifying epigenetic signals into two categories: non-coding gene expression signals and epigenetic memory signals (variable order). In Figure 1, the cell cycle is simply divided into two stages, interphase and telophase. Gene expression occurs during interphase, while cell division and related processes occur in telophase. In interphase, gene expression is controlled by stable epigenetic signals on master regulatory genes. In telophase, the protein environment can alter the epigenetic marks of master regulatory genes [15].
Yoko Suzuki et al. showed that time delays completely alter time-dependent dynamics even for the simplest possible circuits with one or two gene elements by self-regulation and regulation of each other. These elements can cause complex behaviors such as periodic, quasi-periodic, weak chaotic, strong, and intermittent chaotic dynamics. They introduce a special power spectrum-based method for describing and distinguishing these dynamic states. They argued that cancer could cause further delays in the regulation of circadian genes, thus creating a non-periodic dynamic that disrupts the normal process. In fact, the circadian gene regulator circuit can potentially induce chaos in the circuit by delay [16].
DandanYue et al. investigated a fractional order genetic regulatory network and several of its dynamic behaviors. Their numerical simulations showed the effect of fractional order on stability and oscillations. They presented the results of fraction stability by considering degradation rates as a bifurcation parameter to study the dynamics of the system. Their findings confirm that fractional differential and diffusion can better describe gene activities and provide a greater understanding of nonlinear features [27].
Binbin Tao et al. presented a new two-gene regulatory model with a delay parameter using the fractional order system that can better describe the memory and intrinsic properties of genetic regulatory networks. Total delay was selected as the bifurcation parameter and sufficient stability conditions were obtained. In this research, it has been discovered that the network stability interval is inversely proportional to the network fraction order and increasing the order can reduce the critical amount of delay [30]. A very similar work with similar results is also presented by Qingshan Sun et al. [31].
In this paper, for the first time, Figure 1. A simplified form of GRN in the cell cycle based on the VOF system model. In this simple platform, epigenetic factors in the form of signals affect the network of gene interactions. The first is through the expression signals of non-coding RNAs and the second is through epigenetic memory which is equivalent to the variable order signals of the system. The letter P represents a protein and the letter G represents a gene. Dashed lines indicate indirect communication of genes through proteins.
In this paper, for the first time, GRNs are modelled using nonlinear delayed VOF systems via gene expression time series data. Also, the dynamics of a real GRN system that fits on real data are discussed. For this purpose, we have created a new tool called GENAVOS (Gene Network Analysis by Variable Order Systems). In the first step, time series data are entered into the tool. Then by defining the network structure graph, the nonlinear VOF differential equations of the network are generated in the state space automatically. An advantage of this tool is the ability to approximate the variable-order functions of the system over a period of time. This is possible by selecting several methods as function approximation with free parameters, including the use of radial base neural networks, Fourier series, constant function, cosine function, and random function. It is also possible to model networks in which one or more modules do not have time series data such as non-coding RNAs. This is done by approximating the expression of dateless genes with the RBF neural network. In the next step, the numerical solution of the VOF system is performed simultaneously with determining and optimizing the system parameters via minimizing the sum of absolute error [32]. The parameters identification & optimization is done by the imperialist competitive algorithm (ICA), which is an evolutionary optimization algorithm [33]. After implementing the GRN system, GENAVOS offers several tools for analysing the system dynamics. In the continuation of this article, we present the working methods and the results of this tool. GENAVOS is available at https://github.com/hanify/GENAVOS with its user guide and MATLAB codes.

Time Series Dataset
Gene expression time series data of cell cycle process are suitable for analyzing GRNs. We used two sets of data. The first set is related to human primary skin fibroblasts (PSF) throughout the cell cycle, which is a normal cell [34]. This dataset is available under the access number E-TABM-263 at https://www.omicsdi.org/dataset.
The second is the human cell cycle in the Hela cell line, which is a cancer cell [35]. This dataset is available at http://genome-www.stanford.edu/Human-CellCycle/HeLa/. Unfortunately, despite the vast amount of gene expression data, the number of gene expression time series data that have the desired number of time points and describe a biological process such as cell cycle is small. It is so frustrating for the non-cancerous cells that we did not find any suitable data other than primary skin fibroblasts data. Hela cell cycle data have several time series with different number of time points. However, because the primary skin fibroblasts data are related to a time series with 2-h intervals in the cell cycle, we also chose a time series with 2-h intervals for the HeLa cell.

Case Study and Network Structure
Because the selected time series data are relevant to the cell cycle, we focus on one of the infrastructures of gene regulatory networks in the cancer cell cycle. In this section, we consider the structure of the cancer gene network regulated by miR-17-92 cluster during G1/s transition in the mammalian cell cycle [36]. This network is shown in Figure 2. In this graph, T-edges indicate gene inhibition and arrow edges indicate gene activation. This graph is equivalent to adjacency matrix [a ij ] where a ij = 0 means no connection between node i and j, a ij = 1 means node i is activated by node j and a ij = −1 means node i is inhibited by node j. In GENAVOS software, the adjacency matrix of the network graph must be defined.

GRN Modelling by VOF Systems
The variable fractional order derivative is a generalization of the fractional order derivative. Like the fractional order derivative, several definitions are proposed for the fractional variable order derivative, with Riemann-Liouville, Caputo, and Grünwald-Letnikov being the most common definitions. In this article we have used the left fractional VO derivative of Grünwald-Letnikov. According to this definition, we have [37]: To discretize Equation (1), we assume t i = ih, i = 1, 2, . . . , k, so the numerical approximation of the variable order derivative will be as follows.
Now, if a variable order fractional differential equation is as follows [38], to solve it numerically, we can use Equation (5).
To model GRN in state space with VOF system, we can write the following equation for each node of the network (genes). (6) where F + and F − are nonlinear functions that indicate activation or inhibition of the G i gene by its control genes, respectively. θ is the parameter of nonlinear functions F,γ is the degradation parameter, κ is the synthesis rate parameter, and τ is the delay parameter. The j index represents the genes that activate or inhibit the G i gene and is determined from the adjacency matrix [a ij ] of the network graph [39]. Although F functions can be any nonlinear function, experimental studies have shown that they have a sigmoidal form. The two functions known as the sigmoidal form are the Hill Type and Log-Sigmoid functions, which can be selected in GENAVOS software [40]. Therefore, the functions F + and F − can be in the following forms.
Now, to numerically solve Equation (6) using Equation (5) by considering the delay parameter, we first consider the following conditions: Therefore, the final equation for numerical solution (discrete form) will be as follows:

Function Approximation for Variable Order
To solve Equation (6), it is necessary to determine the function of the variable order, α(t). In real systems, α(t) is a time-varying variable (parameter) that affects the system In the GRN system, for example, this is most likely a time-varying epigenetic signal. For phenomena such as biological systems, it is almost impossible to determine α(t) analytically. So, in the optimistic case, α(t) is obtained through measurements as a time series of epigenetic factors such as DNA methylation. However, if the values of α(t) are measured or sampled with the time step of solving equation (h), Equation (6) can be solved numerically and not analytically. However, in fact, the measurement intervals of an experiment are greater than the time step h. So, we have to approximate α(t) from the measured values (time series). GENAVOS has several solutions for determining α(t), including determining α(t) as a constant function, a function with random values, or approximating it by the Fourier series, radial base functions (RBF neural network), and the Cosine function. The best results were obtained with RBF neural networks. As we all know, RBF neural networks can approximate any function with acceptable error [41]. Therefore, α(t) is approximated as follows: where c is the function parameter and the center of the Gaussian kernel, ρ is the coefficient of the Gaussian kernels, and L represents the number of Gaussian functions. The second part of Equation (11) guarantees that 0 ≤ α(t) ≤ 1.

Expression Approximations for Genes That Have No Data
Gene-level regulatory elements include non-coding and coding genes, as seen in Figure 2, which contains the miR-17-92 cluster. Unfortunately, time series data similar to those for coding genes are not available for miRNAs and LNC-RNAs. Here, we have included an initiative in GENAVOS to approximate the expression of genes that lack data. The procedure is very similar to the order function approximation described in Section 2.4. Here, again, an RBF neural network is considered for the expression of the no-data gene (here miR-17-92) according to Part I of Equation (11). We show the parameters of this network with ζ and µ so, we have:

Parameters Identification and Optimization
The purpose of the parameter estimation step is to determine and optimize all system parameters so that the system's solved response (here. predicted values of gene expression) has the least error with real data (here, gene expression time series data). In GENAVOS this is done by the imperialist competitive algorithm (ICA). ICA is a method in the field of evolutionary computing that finds the optimal answers to various optimization problems. In terms of application, this algorithm falls into the category of evolutionary optimization algorithms such as genetic algorithms, particle swarm optimization method, ant colony algorithm and so on. Like all similar algorithms, the ICA forms the initial set of possible answers. These initial answers are known as "country". ICA, using the assimilation, revolution, and imperialistic competition operators, gradually improves these initial answers (countries) and finally provides the appropriate answer to the optimization problem [33].
To use this algorithm, the sum of absolute error (SAE) function is defined as a cost function as follows: where N s is the number of time samples of gene expression data, N g is the number of genes in the network, and T s is the time step of data sampling. The parameters we are trying to optimize are: {κ, γ, θ, c, ρ, µ, ζ}.

Summary of the Proposed Method
Epigenetic signals are modelled by RBF neural networks as shown in Figure 1, although other methods and functions can be used in GENAVOS. These signals are divided into two categories: (1) α(t) as the variable order of the system and the agent that creates time-varying memory, and (2) the expression of non-coding genes. The ICA method, as an iterative algorithm, determines the optimal parameters of the VOF system plus the parameters of the epigenetic signal models. Each time the ICA repeats, the output of the VOF system, which is the time series signals of the gene expression, is compared with the actual data by calculating an SAE error function (or any other appropriate error function). This error is returned to the ICA algorithm as the objective function. The reason for using SAE instead of other error functions, such as MAE or MSE, is its high error resolution because it does not use averaging to reduce its numerical value virtually. Figure 3 shows the workflow of this method.

GENAVOS Analytical Tools for System Dynamics
After implementing the system and identifying its parameters, we can study the dynamics of the GRN system using GENAVOS analytical tools. The first is related to sensitivity to initial conditions. Here, by applying a coefficient (for example, 1.01 or 0.99), the initial conditions are slightly changed to determine whether the system error changes in a limited way or not. In other words, has the system relatively stable response under these initial conditions? Note here that although obtaining the equilibrium point of the system described by Equation (6) is straightforward (obtained by equating the other side of the equation to zero), but examining the stability of the system using the Jacobin matrix characteristic equation is not straightforward due to the existence of a time-varying order. As we know, due to the delay parameter, we need to take the Laplace transform to obtain the characteristic equation of the linearized system [30]. For the fractional order derivative, the Laplace transform is definable as well as the integer order. However, for the variable order derivative, the definition of Laplace transform is only possible if α(t) is explicitly expressed as a function [42]. The second is to draw the bifurcation diagram for changing system parameters. These graphs can show which state variables (here, genes) change with parameter change and how the system exhibits periodic, quasi-periodic, chaotic, and other behaviors. This diagram can show which parameter change derails the GRN system and leads to disease, for example by changing the delay, degradation, or synthesis rate parameter. It is also possible to check the system's temporal response to changing parameters.
The third is the chaos test with the 0-1 algorithm. This test is used to distinguish regular state from chaotic state in dynamic systems. Unlike the Lyapunov maximum exponent test, which uses phase space reconstruction, the 0-1 test uses a time series of data obtained from a dynamic system [43]. The test input is a one-dimensional time series ϕ (n) for n = 1, 2, ... We use ϕ (n) data to construct the following 2D system: p(n + 1) = p(n) + ϕ(n) cos(cn) q(n + 1) = q(n) + ϕ(n)sin(cn) n = 1, 2, . . . , N where c ∈ (0, 2π) is a random number. The bounded trajectory on the q-p plane indicates a regular system, and the trajectory, similar to Brownian motion or random walking, indicates a chaotic system. The (time-averaged) mean square displacement is defined as follows: and its growth rate is as follows: R is the test result that can be determined by two methods of regression and correlation coefficient [43]. In this work, we have used the correlation coefficient method. According to this method, if R is close to zero, it indicates a regular system, while if it is close to one, it indicates a chaotic system.

Results
This section summarizes the results. These results are based on non-cancer cell data to simulate how a normal cell system can turn into cancer. Similar results are given for cancer cells in Supplementary Results. Table 1 shows the initial parameters and initial settings of the system in accordance with what we have already described. ICA settings are listed in Table S1. The first time point of the data is as control data at t = 0. Therefore, we consider these points as the initial conditions of the equations of state space.  Table 2 shows the value of the SAE function in the two modes α(t) = 1 and α(t) approximated by RBF networks. As can be seen, the SAE differs significantly between the two modes. Table 3 shows the effect of the delay parameter on SAE, so the best results are obtained for τ = 0.1 and τ = 0.3 for normal and cancer cells, respectively. Figures 4 and 5 show the results of numerical solution of the normal and cancer cell systems with real data points at 2 h intervals. As shown in Figures 4 and 5, the modelled signals fit the data very well and mimic the behaviour of the data. Figure 6 shows the order function curve for each system state variable. Table 4 shows the values of the obtained parameters for the synthesis rate and degradation. There is a similar table for Hela cells in Table S2.     Figure S3 shows the system response to changing initial conditions. In this figure, the initial conditions are up and down with a coefficient of 1.1 and 0.9. Time delays can cause oscillatory gene expression and provide insights into the dynamics of interactions between genes associated with cancer suppression [6]. Therefore, the delay parameter was selected as a bifurcation diagram parameter. In Table 4, the degradation and self-synthesis rate of miR-17-92 is zero. Because the role of miR-17-92 in the cancer cell cycle is very effective, we chose the degradation rate of miR-17-92 as one of the bifurcation parameters. It has also been shown that a gene with self-inhibiting feedback causes oscillations in the system. Therefore, we chose the MYC synthesis rate as another parameter for the bifurcation diagram. Figure 7 shows the system bifurcation diagram for changing several parameters. The system exhibits complex dynamics similar to chaos or complex oscillations for ranges of parameters.  Figure 7A,B shows the behavior change of the two system state variables (CDC25A, E2F1) for the delay parameter change. Figure 7E shows the evolution of the MYC gene to change the self-synthesis rate parameter. Figure 7F shows the change in behaviour of E2F1 to change the miR-17-92 degradation rate. Figure 7C,D shows the effect of the variable order α(t). To achieve this goal, we have defined the Q parameter as Equation (17):

System Dynamics Evaluation
where α Normal (t) is obtained from modelling normal cell data (Primary Skin Fibroblasts) and α Cancer (t) is obtained from modelling cancer cell data (Hela Cell) by the proposed method. Figure 8 shows the modelled time signals corresponding to the bifurcation diagram of Figure 7. Finally, to check for chaotic dynamics in the model, we used the 0-1 test [43]. The q-p plots of normal and cancer cell systems are shown in Figures S9 and S10. Table 5 shows the results of this test in different system conditions. The q-p plots corresponding to Table 5 are shown in Figure 9.    Table 2 shows that α(t) has created more flexible dynamics in the time-varying mode. It seems very difficult or impossible to fit the system response on the data without a variable-order differential, while Figures 4 and 5 show that the VOF system fits well on the data. Although previous works (such as the mentioned in the introduction) have provided flexible models by delayed fractional-order systems, none of them have implemented the system on real data. In GENAVOS software, we have made it possible to take the fractional order as a fixed number, similar to what have been done in previous works. Figures S1 and S2 show the results of this selection. As it turns out, without variable order, the model has a very weak fitness on real data. Of course, our goal in this work was not to build a model that fits well only on data because many other reverse engineering models such as neural networks and nonlinear regressions can also fit well on data. Rather, our goal was to build a model based on the biological reality of GRNs with the laws of biophysics and biochemical kinetics that are consistent with real data so that it can justify or predict GRN-induced cancer cell phenomena. It is clear that the data may also have noise values, however, several articles have explicitly emphasized the existence of oscillating dynamics and periodic behaviors of gene expression data in the cell cycle, and this has been biologically proven [34,35,44]. Thus, fluctuations in gene expression over the course of the cell cycle are largely related to the intrinsic properties of the system. In this regard, Mahboobeh Ghorbani et al. studied the mathematical characteristics of gene expression dynamics in GRNs using the Hurst exponent [45]. The Hurst exponent H is used as a measure of long-term memory (also called long-range dependence LRD) of time series. The value of H is between 0 and 1, so the value of 0.5 indicates the absence of long-term memory. A value greater than 0.5 indicates a larger degree of long-term memory. They suggested that gene expression is not random, because the Hurst exponent for most of the genes is significantly higher than 0.5. So the gene expression time series cannot be regarded as a random process [45]. This indicates long-range dependency for gene expression signals. The existence of long-term memory for gene expression signals is well accepted due to its slow transcriptional and epigenetic dynamics and cellular memory. To evaluate the fit of the proposed model on this nature of the gene expression time series data, we calculated the Hurst exponent for the modelled gene expression signals given in Table S3. As shown in Table S3, H is significantly greater than 0.5. It has also increased from normal to cancer for all genes except the miR-17-92 cluster. This means an increase in the degree of long-term memory in the cancer cell compared to the normal cell. Since a major factor in creating system memory in the VOFS model is the variable order α(t), we expect the same behavior for the average value of signals α(t). Table S4 shows the average value of α(t) for each state variable. This table shows that the mean value of α(t) decreased from normal to cancer, except in the case of miR-17-92. Therefore, the proposed model fits well with the data and their mathematical nature.

Realistic Prediction of the Role of the miR-17-92 Cluster
The miR-17-92 cluster in the human genome encodes seven mature microRNAs. It has been validated as regulator with both oncogenic and tumor suppressive properties [46]. Many studies show that the miR-17-92 cluster is overexpressed in cancer. This is specifically stated in a review article by Xinju Zhang et al. [47]. As shown in Figure 4, this module is under-expressed in normal cells. Instead, it can be seen in Figure 5 that its expression has increased significantly and has fluctuating states during the cell cycle. Table 4 shows that the self-synthesis and degradation parameters for the miR-17-92 cluster are zero. This means that this module (resulting from the activity of seven miRNAs in the cluster) is inactive in a normal cell. While the values of the mentioned parameters in Table S2 show that this module is activated in the cancer cell. This is consistent with biological studies [46,47]. Another interpretation of miR-17-92 is its different Long-term Memory behaviour in normal and cancer cells according to Tables S3 and S4. A justification for reducing the degree of long-term memory of this module in cancer cells relative to normal cells (or the faster dynamics in a cancer cell than in a normal cell according to Table S4) can be as follows: The presence of long-term memory is most likely due to epigenetic memory factors, which in our model are α(t) signals. The miR-17-92 module is also an epigenetic factor. Therefore, it is affected by epigenetic memory differently from the coding genes. Therefore, our model approximates the function of miR-17-92 as an epigenetic agent in accordance with its biological reality. Our model could be an approach to simulate and synthesize gene expression time series for non-coding RNAs.

The Role of System Parameters in GRN Dynamics
The role of three parameters in GRNs is key: delay, synthesis rate and degradation. Regulation of gene expression involves the synthesis of mRNA and protein through transcription and translation and degradation of molecules [48]. In other words, genes in a gene regulatory network influence the rate of synthesis and degradation of each other. Translational and transcriptional delays, as well as other environmental and epigenetic factors, also affect the rate of synthesis and degradation [49]. Therefore, in general, these three parameters can play two basic roles in regulation of gene expression: (A) they control gene expression levels (steady state mRNA concentration). In other words, they control the bias level of gene expression signals. (B) They control the fluctuations of gene expression signals on the bias line. It is noteworthy that with the change of these parameters, the expression of the three genes, MYC, E2F1, and CDC25A, undergo serious fluctuations. They are known to be oncogenes and potential suppressors of tumours, and fluctuations in their expression are crucial in cancer. Table 3 confirms that the delay parameter in cancer cell GRNs is larger than in normal cells. This confirms studies that have shown that cancer increases the delay parameter [16]. The increase in the delay parameter in cancer was not unexpected due to the high complexity of the cancer network due to increased interactions with intermediate modules such as proteins.
As we expected, Figure 7 shows that the system undergoes a fundamental change by changing its parameters. The E2F1 gene undergoes complex oscillations from the τ ≈ 5.31 onwards, while returning to regular behavior between τ ≈ [8. 37, 9.27]. The oscillation starts again from 9.27. A similar behavior is seen for CDC25A. Oscillation start at around 2.7 and become almost regular at ≈5.59 and then oscillate again. A definite opinion as to whether the system of these oscillations is chaotic or quasi-periodic is not straightforward by simply using a bifurcation diagram without analytical methods. That's why we used the 0-1 test. Bifurcation diagram of MYC is similar to the bifurcation diagram of twodimensional maps (such as logistic map). This is due to the presence of a negative feedback loop of the MYC gene on itself, which behaves like a recursive map in the simulation. The E2F1 gene oscillates by changing the degradation parameter of miR-17-92 from about 2.51, due to the presence of an inhibitory feedback loop between them.
In Figure 8, the time signals confirm the behaviors of the system in a range of parameters where it is irregular corresponding to Figure 7. A comparison of Figure 9 and Figure S9 clearly shows that changing the parameters corresponding to Figure 8 makes the behaviors of the trajectories irregular and resembling a random state. Table 5 also confirms this. So, the signals in Figure 8 are somewhat chaotic. However, given their R-values in Table 5, this is probably weak chaos. Table 5 also has an important result: R is larger for the cancerous condition than normal. In other words, the weak chaotic property of the system in the cancerous state is more than the normal state.

Control of GRN Dynamics by Possible Epigenetic Signals
We find three important interpretations of Figure 6. The first is the existence of very slow dynamics in regulatory networks. We calculated the average value of α(t) for each state variable whose results are shown in Table S4. We then obtained the sum of these averages for both cancer and normal cells. (In the case of integer order systems, this sum represents the order of the system). Its value was 1.88 and 1.33, respectively, and in general it seems that the dynamics in the cancerous state has slowed down. If we consider α(t) as a system input (and not a system variable [42]) that is affected by epigenetic memory, this means slowing down epigenetic mechanisms in cancer [25]. The second is the frequent zeroing of the α(t) function over time, which evokes the type of on-off switching mentioned in the work of Richard P. Halley-Stott et al. [26]. It may also support studies that claim that gene expression is inherently discrete [50]. Another conclusion from α(t) is their role in controlling the temporal pattern of gene expression as possible epigenetics signals. In the biological context, potential epigenetic signals are signals generated by epigenetic memory mechanisms such as DNA methylation and chromatin remodelling, as described in Figure 1. For this conclusion, we extracted the extremum points of signals α(t) and gene expression signals and then plotted them together in two curves, as shown in Figure  10. We found that a minimum/maximum point of α(t) results in a maximum/minimum point for gene expression signals after a delay (the maximum occurrence in one leads to the minimization of the other, and vice versa). The amount of delay between these points is not constant and can be greater than the total system delay, i.e., τ. It varies between 0.1 and 0.5 for both normal and cancer cells. An important reason for this delay being larger than τ is that epigenetic dynamics are slower than transcriptional dynamics. Now, a question arises: Is it possible to simulate the behavior of a cancer cell by changing the simulated epigenetics of the normal cell system? To answer this, we first change the delay parameter of the normal cell system from 0.1 to 0.3. Then in Equation (17), we take Q equal to 1. When Q = 1, α new (t) will be equal to α Cancer (t). In the last step, we change the synthesis rate and degradation of miR-17-92 from zero to 1.86 and 0.9, respectively, according to Table 4 and Table S2. Figure 11 shows the time response of the system to substituting α Normal (t) with α Cancer (t). We found that the normal cell system shows exactly the same behavior as a cancer cell by replacing α(t) with α Cancer (t). In other words, the time patterns of gene expression in normal cells are similar to their time patterns in cancer cells under the influence of signals of variable-order functions, but the levels of gene expression are different. This indicates that the VOF-GRN system is strongly influenced by the variable-order function. The order functions are not parameters that relate to just its own state variables, but system inputs that affect the system as a whole. In biological terms, transcriptional kinetic parameters (such as the rate of synthesis and degradation affected by the expression of other genes and proteins) probably control the expression level of a gene in GRNs, while possible epigenetic signals control its temporal pattern. Figure 11. The normal cell mimics the temporal pattern of gene expression in cancer cells by substituting α Normal (t) with α Cancer (t) as the epigenetic factor.

Conclusions
Recently, the theory and application of variable-order differential equations has expanded. These equations are very useful for describing flexible and reality-based systems. In this paper, a new tool called GENAVOS was introduced to model gene regulatory networks. Our results show that the use of delayed nonlinear VOF systems has very flexible dynamics that correspond to real biology. The VOF system for GRN is so strongly influenced by the time-varying order function as possible epigenetics signals that we conclude that variable-order functions as potential epigenetic signals probably control the temporal pattern of gene expression, while genetic factors control the level of gene expression, and the two factors are constantly interacting with each other. It is suggested that the variable order functions obtained as possible epigenetic signals in this study be compared or replaced with epigenetic time series data such as DNA methylation.
Supplementary Materials: The following are available online at https://www.mdpi.com/2073-899 4/13/2/295/s1, Figure S1: Modeled gene expression signals for PSK cell with α(t) = 0.8, Figure S2: Modeled gene expression signals for Hela cell with α(t) = 1, Figure S3: Normal cell system response to changing initial conditions. In this figure, the initial conditions are up and down with a coefficient of 1.1 and 0.9, Figure S4: Variable-order functions of cancer cell system approximated by RBF neural network, Figure S5: Phase plane diagrams of normal cell system, Figure S6: Phase plane diagrams of cancer cell system, Figure S7: Bifurcation diagram for changing MYC Self-Synthesis Rate in cancer cell, Figure S8: Cancer cell system response to MYC Self-Synthesis Rate κ 1 =50, Figure S9: Q-P plots for normal cell, Figure S10: Q-P plots for cancer cell, Table S1: Imperialist Competitive Algorithm Parameters, Table S2: Values of the obtained parameters of normal cell system, Table S3: Hurst Exponents of modelled gene expression signals (Time series), Table S4: Average value of variable orders.