A control theoretic three timescale model for analyzing energy management in mammalian cancer cells

Graphical abstract


Introduction
Cellular decision making and responses are orchestrated by a set of complex biochemical pathways/networks. Broadly speaking, biochemical pathways/networks can be categorized as metabolic pathways, gene regulatory networks (GRNs), and signaling pathways. A metabolic pathway is a coherent set of biochemical reactions catalyzed by a number of enzymes. It helps a living organism to transform an initial (source) compound into a final (target) compound and energy. On the other hand, fundamental information processing and control mechanisms in a cell are performed by GRNs. Regulatory genes code for proteins that activate or inhibit the expression of other genes. Thus, a complex web of interactions, called a GRN, in terms of activation and inhibition of genes, is formed. Moreover, signaling pathways contain a series of specific actions in a cell in which a signal is passed from the environment and accordingly the cell responds. These pathways are of diverse nature and are interacting with one another. They also form a hierarchy. For example, an entire organism can be thought of as a huge network of interacting organs, each of which, in turn, is a network of interacting tissues. A tissue is a collection of interacting cells performing similar functions. A cell may be considered as a huge network of interacting components to constitute aforesaid biochemical networks. Moreover, the influence of environment and other factors on enzymes or gene regulation needs to be considered to make a study on how an organism is responsive to environmental changes. Therefore, it is necessary to focus not only on individual processes or pathways but also on their integration. among metabolic, signaling and gene regulatory networks related to central carbon metabolism in a mammalian cell. Attempts to elucidate such huge biochemical networks on the structural/ functional basis face the problem of combinatorial explosion. There exist several investigations on modeling each of these pathways individually [1][2][3][4][5][6][7][8][9]. However, the study on the integration of gene regulation, metabolism, and signaling events (pathways) is not much.
Among several approaches one of the most common is to explore metabolic pathway is Flux Balance Analysis (FBA) [10,11]. FBA can analyze only steady-state response of a system rather transient response. However, Metabolic Control Analysis [12] overcomes the limitation of FBA. Yet, it does not provide any supervisory controller that can regulate the enzyme/metabolite concentration to meet some particular requirements of a cell. Recently, we have developed a control theoretic approach to solve this problem [13][14][15]. Some investigations on crosstalk mechanism [16] and identifying properties [17,18] of signaling networks have also been developed. Since the activation of regulatory mechanism depends on corresponding gene expression level thus any GRN can be modeled mathematically by coupled ordinary differential equations (ODEs), Boolean network or Bayesian network [19].
Response time of the interconnected complex subsystems may vary from each other. In other words, they can be treated as complex subsystems of different timescales. Biochemical pathways, such as metabolic, signaling and regulatory networks show this kind of phenomena. Timescale differences among metabolic, signaling and gene regulatory networks make their integration a daunting task. Several researchers have tried to model the dynamic behavior of a system by stoichiometric reconstruction of integrated biochemical pathways [20]. Here, the authors have considered two different timescales. The metabolic and signaling networks have been considered in fast timescale, whereas GRN has been considered in slow timescale. However, the metabolic, signaling and gene regulatory networks participate in three different timescales. Previous investigation [21] show that proteins in Escherichia coli take seconds to minutes to express. It will be longer in the case of mammalian systems. On the other hand, metabolite concentrations can respond in seconds to microseconds [22]. Gene regulatory events take minutes to hours [20]. Thus, integration of these biochemical pathways becomes a three timescale problem, where metabolic, signaling and gene regulatory networks are ultrafast, fast and slow respectively.
In this study, we develop a three timescale multiple input and multiple output (MIMO) model based on support vector regression (SVR) and genetic algorithm based controller to simulate the dynamic behavior of integrated signaling, metabolic and gene reg-ulatory networks specifically responsible for mammalian carbon metabolism in normal and cancer cells. Here, the metabolic, signaling and GRN have been mathematically modeled by sets of ODEs with proper three timescale selection, feedback, allosteric effects and perturbation. In this context, we have collected all pathway information under consideration from KEGG database [23] and literature. Finally, we have developed a genetic algorithm (GA) based controller to drive the change of concentrations/expression levels of certain metabolites/proteins/genes with respect to time in a desired fashion. Thus, the proposed model can also predict the possible effects of certain drug targets of specific diseases through GA controller.
The central carbon metabolic (CCM) network comprises glycolysis, tricarboxylic acid (TCA) cycle, pentose phosphate pathway (PPP) and free fatty acid (FFA) metabolism. Glycolysis consumes glucose to produce pyruvate and energy in the form adenosine triphospate (ATP), although most of the energy is generated from TCA cycle taking pyruvate as input. Here, PPP is involved in macromolecular synthesis including reduced nicotinamide adenine dinucleotide phosphate (NADPH) and ribose 5P. FFA is formed by acetyl CoA of TCA cycle. However, a mammalian tumor cell slows down oxidative phosphorylation by inhibiting TCA cycle activities and produces abnormal amount of lactate. This effect is called ''Warburg effect" [24,25] which may lead to cancer [26]. In this context, we have analyzed the energy and cell proliferation management of CCM pathway along with corresponding signaling and gene regulatory networks in mammalian cancer cells. Here, we have considered the interactions among different enzymes/proteins, transcription factors and genes associated with central carbon metabolism to capture the cellular dynamics during normal and malignant conditions in mammalian cells. These detailed interactions can be found in Supplementary Tables S1-S8 (abbreviations of different molecules can be found in Supplementary Table S10). Fig. 1 depicts the integrated metabolic, signaling and gene regulatory networks associated with central carbon metabolism. We also have analyzed the possible effects of six drug targets, such as deactivation of pyruvate kinase, glucose-6-phosphate dehydrogenase, transketolase, ribose 5P isomerase, glucose-6-phosphate isomerase and finally activation of pyruvate kinase, on mammalian malignant cells.
We have also compared the proposed model with the following theoretical models. In this context, FBA based model [66] has depicted a similar set of malfunctioning enzymes as identified by the proposed model in cancer cells for its growth. However, FBA based model does not involve transient analysis of the behavior of molecules in cancer pathways. Previous differential equation and optimization-based models [67][68][69] have also depicted similar mutation as depicted by the proposed model in cancer cells. However, these models do not consider the three timescale nature of the integrated CCM pathways. Besides, they have not explored the outcome of probable drug targets to control energy metabolism in cancer cells.
Here the proposed model is able to predict that deactivation of glucose-6-phosphate dehydrogenase and ribose 5P isomerase may slow down growth and proliferation of cancer cells. However, glucose-6-phosphate dehydrogenase deactivation may not be able to reduce fermentation and energy supply in malignant cells. On the other hand, though ribose 5P isomerase deactivation may inhibit cell fermentation, it may not get success to stop energy supply in mammalian cancer cells. In this context, deactivation of transketolase and glucose-6-phosphate isomerase may be considered as potential drug targets to resist cell growth, fermentation and proliferation as well as energy supply in human cancer cells. Activation of pyruvate kinase (M2 isoform) may also reduce cancer progression in terms of cell growth, proliferation and fermentation. Besides, it may fail to decrease ATP production in tumor cells. However, deactivation of pyruvate kinase may be a poor choice as a drug target for cancer therapy.
In summary, the proposed novel model can tackle three timescale nature of any integrated biochemical pathway comprising metabolic, signaling and gene regulatory networks. Subsequently, it can capture the nonlinear transient dynamics of the integrated network of both normal and perturbed human cells. This has been validated appropriately by some other methods/results [70][71][72][66][67][68][69]. Moreover, using proposed GA controller, effects of drug targets on diseased cells have been explored and analyzed. Thus, it can help to find out a novel therapeutic goal for complex diseases, such as cancer and type 2 diabetes.
Number of chromosomes and generations respectively.

Method
Here, we describe the proposed methodology for developing the model that will mimic the behavior of an integrated pathway system. The methodology involves six steps. At first, timescales of metabolic, signaling and gene regulatory networks have been selected, and the basic dynamics of the integrated pathway system has been described. Then, we have formulated the state equations for the integrated pathway. In the third step, an appropriate MIMO plant has been developed by tuning its parameters. Subsequently, training dataset has been generated from the appropriately perturbed MIMO plant/system. Thereafter, SVR model is developed to approximate the MIMO plant. Finally, GA controller is applied on the plant as a model predictive controller to capture nonlinear transient cellular dynamics. Fig. 2 depicts the flowchart of the entire methodology. Here we have considered central carbon metabolism related metabolic (glycolysis, TCA cycle and pentose phosphate pathways as well as FFA synthesis/consumption), signaling and gene regulatory networks throughout the entire methodology. All the mathematical symbols used in this article are defined in Table 1. However, each of these symbols has again been described whenever it appears first.

Selection of three timescales and the basic dynamics of an integrated pathway system
In an integrated biochemical pathway system, involving metabolic, signaling and gene regulatory networks, the rate of a metabolic reaction is controlled by an enzyme/protein, which is originated from an expressed gene participating in a gene regulatory network. The expression of this gene is in turn regulated by one or more transcription factors (proteins) which are expressed/ activated by signaling molecules in a cascade of signal transduction pathway. Here fð:Þ; gð:Þ and hð:Þ represent activities corresponding to gene regulatory, signaling and metabolic networks. The time (t) derivatives of x; y and z are represented by _ x; _ y and _ z respectively. State variable x is slow, whereas y and z are fast and ultrafast state variables respectively. The symbols x 1 and x 2 correspond to slow to fast and fast to ultrafast timescale ratios. Here, x 1 and x 2 are small positive quantities. According to previous investigation [21], proteins in Escherichia coli take minutes to express. It will be longer in the case of mammalian systems. However, metabolite concentrations can alter in seconds [22]. On the other hand, gene regulatory events take hours [20]. Based on these assumptions, we have considered x 1 ¼ 1=60 and x 2 ¼ 1=60, and thereby the stretched timescales as s 1 = t=x 1 and s 2 = s 1 =x 2 = t=ðx 1 x 2 Þ. Thus the timescales for gene regulatory, signaling and metabolic networks are t; s 1 and s 2 respectively.
We now obtain the fast (ĝð:Þ) and ultrafast (ĥð:Þ) subsystems corresponding to signaling and metabolic networks from the top down timescale decomposition [73,74]  Again, the fast and ultrafast subsystems corresponding to signaling and metabolic networks can be treated as a two timescale system. The ultrafast subsystem corresponding to metabolic network (ĥð:Þ) is given by ðx; y; z; uÞ ð 14Þ Here we can determine the value of u, in terms of x; y and z, considering the equilibrium state of the system.

Formulation of state equations of integrated biochemical pathways
Here we have shown how state equations for gene regulatory, signaling and metabolic networks are formed. Let a gene x i be expressed by s 0 (1 6 s 0 < s) transcription factors y j 0 (j 0 ¼ 1; . . . ; s 0 ) included in a signaling pathway. Thus, _ x i can be written as The terms e i and b i represent the expression and basal production (of mRNA) rates [76] of x i respectively, whereas d i denotes the decay rate. The expression rate e i of a particular gene x i depends on the binding of many required transcription factors. Absence of any of those transcription factors, the expression of that particular gene is driven by its basal production and decay rate only. Now, e i can be defined as where K ðgÞ ij 0 is the binding rate constant for j 0th transcription factor required to express the gene x i . Let us consider m 0 (1 6 m 0 < m) metabolites z k and c 0 (1 6 c 0 < c) external inputs u l acting as the cofactors [77] to activate an i th gene and Q s 0 j 0 ¼1 y j 0 > 0. As a result, the expression rate e i of the gene x i is enhanced by a multiplication factor. Thus, _ x i can be rewritten as Here F ðgÞ ik and F ðgÞ il represent the constants corresponding to z k and u l to take care of the relative strength of binding with i th gene [13,78,14]. We have considered binding constants to be in ½0; 1Þ. Zero binding constant signifies no binding corresponding to a molecule. Whereas, higher value of a binding constant indicates stronger binding corresponding to the molecule. On the other hand, if m 0 metabolites z k ; s 00 (1 6 s 00 < s) signaling molecules/transcription factors y j 00 (j 00j 0 ; j 00 ¼ 1; . . . ; s 00 ) and c 0 external inputs u l slow down the activation of the i th gene and Q s 0 1 y j 0 > 0, the expression rate e i is decreased by a fraction. Thus, Eq. (17) can be modified as The term F ðgÞ ij 00 represent the binding constant corresponding to y j 00 with i th gene.
A signal transduction pathway can be defined by its interaction matrix N ðsÞ which contains the information about interactions among the signaling molecules including transcription factors. If there are q interactions with initial rate vector r ðsÞ (in timescale s 1 ) of dimension q involving s signaling molecules including transcription factors, the order of N ðsÞ becomes s Â q. It may be mentioned here that r ðsÞ is analogous to a flux vector in a metabolic pathway. The ðj; rÞ th element of N ðsÞ becomes À1 if j th signaling molecule participates in r th interaction. If j th signaling molecule  is produced/activated by r th interaction, the ðj; rÞ th element of N ðsÞ becomes þ1. Let a signaling molecule y j be activated by certain signaling molecules y j 0 (j 0j) in a signaling interaction with initial rate r ðsÞ r ; r th (r ¼ 1; . . . ; q) element of signaling interaction rate vector r ðsÞ . Consequently, the initial interaction rate r ðsÞ r depends on the binding of such signaling molecules activating y j . If some metabolites z k , signaling molecules y j 00 (j 00j 0 ) and external inputs u l slow down the activation of y j , a certain fraction can be multiplied with the initial interaction rate. Thus, we can write using modified mass action law [78] as jj 00 y j 00 Þ where r ¼ rðjÞ is the index for the reaction through which y j is activated. The term K ðsÞ jj 0 represents the strength of binding of y j 0 with y j , whereas F ðsÞ jk ; F ðsÞ jl and F ðsÞ jj 00 represent binding constants corresponding to z k ; u l and y j 00 for the interaction activating y j . Again, if some metabolites z k and external inputs u l accelerate the activation of y j , the initial interaction rate is increased by a multiplication factor. Thus, Eq. (19) becomes Thus, based on Eq. (13), the differential equation representing the dynamics of the signaling network can be rewritten as Similar to the signal transduction pathway, a specific stoichiometry is also associated with a metabolic pathway. Thus, we can consider a stoichiometric matrix N ðmÞ of order m Â n for the metabolic network under consideration, where m is number of metabolites participating in n metabolic reactions in the network with metabolic flux vector r ðmÞ (in timescale s 2 ) of dimension n. The ðk; qÞ th element of N ðmÞ becomes þ1 if k th metabolite is produced through q th reaction. On the other hand, if k th metabolite is consumed in q th reaction, the ðk; qÞ th element of N ðmÞ becomes À1.
Let us consider a metabolite z k be produced by m 0 (1 6 m 0 < m) substrates z k 0 (k 0k; k 0 ¼ 1; . . . ; m 0 ) in a metabolic reaction with initial rate r ðmÞ q ; q th (q ¼ 1; . . . ; n) element of r ðmÞ . This reaction is catalyzed by an enzyme/protein E q (representing the expression level) produced from an expressed gene. Furthermore, let us assume that m 00 (1 6 m 00 < m) metabolites z k 00 (k 00k 0 ; k 00 ¼ 1; . . . ; m 00 ) and c 0 external inputs u l accelerate (noncompetitively or allosterically) the production of z k . Thus, the initial reaction rate r ðmÞ q is increased by a multiplication factor. Then, we can write using modified Michaelis Menten kinetic equation [13,78,14] as where q ¼ qðkÞ (could be one-to-many mapping) is the index for the reaction through which z k is produced. metabolite, represent feedback constants corresponding to z k 00 and u l . Again, if m 00 metabolites z k 00 and c 0 external inputs u l slow down (noncompetitively or allosterically) the production of z k , the initial reaction rate r ðmÞ q is decreased by a fraction. Thus, we can modify Eq. (22) as Now, based on Eq. 14, the differential equation representing the dynamics of the metabolic network can be rewritten as

Design of a MIMO plant
The initial concentrations/expressions of metabolites, signaling molecules and genes have been considered in (0, 1) randomly (Supplementary Table S9). Moreover, the kinetic parameter values have also been initialized randomly within the same interval (Supplementary Tables S12-S15). Based on these initial values, we have simulated the proposed three timescale state space model of an integrated biochemical pathway by solving ODEs in Eqs. 17 (or 18), 21 and 24 with proper timescale selection, and solved for u based on x; y and z values under equilibrium condition. If the model fails to mimic the known behavioral pattern of the integrated biochemical pathway under consideration, we have altered slightly the initial values of kinetic parameters in (0, 1) with the help of previous knowledgebase by trial and error until satisfactory known behavior is captured. Once the proposed model has mimicked the known nonlinear dynamics of the integrated biochemical pathway (as described in SubSection 3.1), we have considered it as normal MIMO plant with the selected values of initial concentrations/expressions and kinetic parameters. This MIMO plant mimics the normal behavior of the concerned integrated pathway. In this context, we have checked whether trial and error based parameter estimation resembles other previous method. For this purpose, we have compared the values of 35 kinetic parameters considered here with those estimated by a method based on hybrid extended Kalman filtering [70]. Fig. 3 depicts that the parameter values in (0, 1) considered in this article are quite similar with those estimated by hybrid extended Kalman filtering based method [70]. To estimate these 35 kinetic parameter values in (0, 1) using hybrid extended Kalman filtering, time dependent true values in (0, 1) of some observed states need to be provided. Consequently, we have used capillary electrophoresis mass spectrom-etry (CE-MS) measured concentration values (normalized in (0, 1)) of eight metabolites in hypoxia-induced metabolic alterations in human erythrocytes [71] as observed states. Here we have externally applied appropriate signal in (0, 1) for corresponding enzyme activities related to hypoxia in human erythrocytes (as described in SubSection 3.1). During the estimation of these 35 parameter values using hybrid extended Kalman filtering, we have noticed that estimated states have been overlapped with observed states as depicted in Fig. 4. However, we could not verify the other parameter values using this way due to lack of required observed states. Now, we perturb the normal MIMO plant by knocking out genes (in this study, genes producing the enzymes pyruvate dehydrogenase and pyruvate carboxylase in carbon metabolism) to capture the nonlinear dynamics associated with a particular biological phenomenon (in this study, ''Warburg effect" [24,25]). As a result, the normal plant becomes the desired perturbed MIMO plant under consideration.

Generation of training dataset
So far we have developed the desired perturbed MIMO plant with consideration of appropriate three timescales. It does not reflect the mutated regulations of a mammalian cancer cell. However, it contains the perturbations associated with ''Warburg effect". In order to capture the altered dynamics of a cancer cell, we need a model predictive controller to be applied on an approximate model corresponding to the MIMO plant under consideration [79]. We are going to develop such an approximate model using support vector regression. Thus, the past input and output data simulated from the MIMO plant under investigation are collected to serve as a training dataset for the SVR-based approximation. Here we have found that consideration of s (=19) number of past inputs and corresponding outputs along with current input (at time h), from the simulated data as an input sample of SVR, and current out-  has been attached to v, as its supervised signal, to form one sample to be included in a dataset D. If there are / such samples in D, let us call it as D / . Here we have considered / ¼ 60000. We have selected 65% data from D / as training data and remaining 35% data as test data to design the SVR model to be discussed in the following section. During training of the SVR model, we aim to utilize a minimum number of a training pattern in order to improve generalization performance. To assess the best possible ratio between train and test data, we have observed the performance with different ratios. The accuracy of such cases is given in Table 2. Here, it can be concluded that the train to test ratio as 65:35 may be the best trade-off between the number of patterns used for training and the test accuracy. Thus we have considered train to test data ratio accordingly.

Design of SVR model
From control theoretic point of view, it is difficult to apply online optimization on actual nonlinear MIMO plant with small sampling periods [80,81]. Moreover, due to large number of vari-ables involved in an optimization problem, and requirement of high sampling rate, it may not be possible to manage online optimization. However, it will be easier if we can provide a model that suitably approximates the nonlinear dynamics of actual MIMO plant. We also think that such implementation may be relevant for future study to determine drug dosages in real time. Even this kind of modelling will be effective in future when actual inputoutput mapping of the MIMO plant is unknown, whereas the inputs and corresponding outputs with time are known. The training dataset D / helps in developing the approximate MIMO SVR  [55] in the case of colon cancer as compared to the normal tissue samples, (D) Similar difference in the expression levels of PI3K in cancer and normal cells have been captured by the proposed computational model, (E) mTOR expression is higher in colon cancer compared to normal one as per the in vivo/ in vitro experiment performed by Zhang et al. [55], and (F) Similar altered behavior of mTOR can be found in our computational results. model to mimic the nonlinear dynamics of the MIMO plant/system under consideration. In this context, it should be mentioned that according to some recent studies [82][83][84], SVR has better performance accuracy than an artificial neural network (ANN) (particularly, multi-layer perceptron (MLP)) and genetic programming (GP). In general, the computational complexity can be evaluated depending mainly on training time. SVR depends on solving a quadratic problem, where a vector and bias need to be calculated. Even during training, only support vectors are selected regarding the kernel. Besides, SVR considers that the number of support vectors is governed by very limited number of training samples. In addition, most of the kernels compute simple dot product. Even SVR may act better than ANN as per the consistency with physical behavior [85]. Moreover, it has been claimed that the SVR technique not only achieves high consistency but also shows a greater degree of generalization to unseen test data compared to some other methods [86]. Besides, a less number of parameters involved in the training phase results in less computation time than that of an ANN [87]. Now, we are going to discuss how a MIMO SVR model can be designed using support vector regression [88,89].
An approximate MIMO SVR model has been developed by combining multiple input single output (MISO) SVR models for each l th Based on Vapnik's e-insensitive loss function, the convex optimization problem corresponding to Eq. (25) can be formulated as [90] min wl;b l ;u;u Ã P l ¼ subject to: c l À w T l dðv a Þ À b l 6 e l þ u a ; for a ¼ 1; Á Á Á ; / w l ; dðv a Þ þ b l À c l 6 e l þ u Ã a ; for a ¼ 1; Á Á Á ; / u a ; u Ã a P 0; for a ¼ 1; Á Á Á ; / Based on the notion of Lagrange multipliers and Karush Kuhn Tucker (KKT) condition [91], the dual corresponding to Eq. (26), in the form of a quadratic programming problem with incorporation of kernel function, can be written as subject to:  [56] has depicted that AKT expression is significantly higher in pancreatic cancer (91 cases) than that in normal pancreas (51 cases). Besides, Roy et al. [57] has shown the deferentially expressed AKT isoforms in normal and malignant oral tissues through immunohistochemical analysis of the human samples, (B) Similar altered behavior of AKT has been shown by the computational results generated from the proposed model.
Here the term jðv a ; v b Þ denotes the kernel function. Besides, Þdðv b Þ ¼ 0 must hold for optimality. Thus, for an a th input, the support vector kernel expansion of Eq. (25) can be obtained fromĉ The summation has been carried out over all the support vec-  Here, we have experimentally chosen radial basis kernel function with kernel coefficient as 10, regularization parameter as 100 and e l =0.0025 to achieve maximum test accuracy.

Design of a controller using a genetic algorithm (GA)
In this study, we have applied a GA [89] on the SVR model as a model predictive controller to control some outputs of the MIMO plant at specific values. There are many reasons behind choosing GA as an optimization technique. Firstly, a genetic algorithm (GA) [92,93] is a stochastic process. Secondly, it is a vigorous search technique, which does not require any information about the structure of the function to be optimized. Such a situation is very common to address a biological system, particularly during mutation. Thirdly, GA is very efficient in handling highly complicated non-linear problems, such as an integrated biochemical pathway system. Besides, due to its inherent parallelism, GA can easily be implemented in a distributed environment. Such a characteristic is very much helpful in handling a large number of parameters/ variables in a biological system. Fourthly, GA may be able to avoid being trapped in a locally optimal solution, unlike traditional methods. Besides, GA can efficiently be applied to different large scale real-world problems with the requirement of multi-objective optimization, such as the particular problem addressed in this article. Such a GA based control mechanism has been depicted in Fig. 2  (Step 6). Here the error between SVR generated predicted output and (user-specified) reference output drives the GA based controller to produce a controlled input based on some cost functions and constraints. The integrated biochemical network (i.e., the MIMO plant) accepts the controlled input to compute actual output. The GA controller has been developed using two algorithms (Supplementary Algorithms S1 and S2). Supplementary Algorithm S1 deals with the working principles of a GA controller. On the other hand, the Supplementary Algorithm S2 drives some outputs of the MIMO plant to target values using Supplementary Algorithm S1 iteratively. Now, we are going to discuss Supplementary Algorithms S1 and S2 in brief. Here, Fig. 5 depicts the flowchart indicating the role and sequence of two Algorithms S1 and S2 for designing GA-based controller.
Supplementary Algorithm S1 shows how proposed GA controller works. Here the GA controller accepts hx; y; z; ui ðinÞ current and v current as input. At first, it generates C chromosomes (set of possible solutions) randomly, forming a population. However, it creates chromosomes around hx; y; z; ui ðinÞ current after getting an optimal solution. For each generation, the fitness value of each chromosome has been computed. Thereafter, the chromosomes are modified by selection, crossover and mutation. In this context, a set of chromosomes is selected according to fitness values greater than a weight value. Subsequently, the crossover is performed by swapping a part of chromosome with corresponding part of another chromosome according to a crossover index. In this context, it should be mentioned that as the size of the output vector is 107 (=ðp þ s þ mÞ as men-tioned earlier), we have considered crossover index as 53. This index is nearly half of the size of the output vector. It is followed by mutation to modify chromosomes with mutation probability 0.7. Here, we have used uniform crossover [94,95], i.e., multipoint crossover, to reduce the mutation bias. These steps are carried out for G generations until the fittest chromosome, whose fitness value is greater than a predefined threshold value, has been evolved as an optimal solution hx; y; z; ui ðinÞ new . Supplementary Algorithm S2 depicts how some outputs of the MIMO plant are controlled to get desired responses with repetitive use of GA controller (Supplementary Algorithm S1). The  [105] has shown that (A) LDH has significantly been up regulated in Dppa4 overexpressed cancer cell, (C) Expression of hexokinase (HK) mRNA levels have been found higher in cancer [106] compared to normal liver tissue, (E) Higher glut1 mRNA levels has been found as 92-fold higher in Meta specimens [106] compared to normal liver tissue. Similar altered behaviors of (B) LDH, (D)  Finally, it should be mentioned that the average value of G is 10 for convergence of Algorithm S1, whereas 30 iterations on average are needed for Algorithm S2 to converge.

Results and discussion
This section has been divided into three subsections. Firstly, we are going to discuss about MIMO plant validation through analysis  [107] in B7-H3 knockdown cells grown in normoxic or perturbed conditions (hypoxia) for 24 h, (B) Li et al. [105] has shown that glucose uptake increases significantly in Dppa4 overexpressed cancer cell, and (C) Our computational results depict similar behavior of glucose uptake in cancer cells compared to that in normal and perturbed ones. of both normal and perturbed behavior of the integrated biochemical pathways under consideration. Secondly, the altered behavior of the integrated biochemical pathways in mammalian cancer cells compared to normal ones will be discussed. Finally, the third subsection deals with the description of the effects of six possible drug targets in mammalian cancer cells.

MIMO plant validation and comparison with other mathematical prediction
Here, we have monitored the concentrations of key molecules during simulation of the proposed three timescale state space model for normal as well as perturbed integrated biochemical pathways related to carbon metabolism. In normal scenario (Fig. S2 in Supplementary material), both glucose consumption and production increase. Glucose is broken down into glucose-6P (G6P). Subsequently, the energy in the form of ATP is consumed. Thus, increasing glucose consumption signifies enhancement of G6P production. The production of ribose 5P and NADPH increases by utilizing higher amount of G6P through PPP. On the other hand, slowing down glyceraldehyde-3P production indicates higher flux through later phase of glycolysis by consuming higher amount of glyceraldehyde-3P. Higher concentration of phosphoenolpyruvate (PEP) supports our claim. Similarly, slowing down the production of pyruvate and acetyl CoA indicates higher flux through TCA cycle and fatty acid synthesis. Enhanced concentrations of oxaloacetate (OAA), citrate and FFA are in conformity with our observation. ATP production increases with higher flux through TCA cycle and later phase of glycolysis. However, after a while ATP decreases due to its higher consumption at early phase of glycolysis. Previous investigations [11,96,97] validate the aforesaid behavior of the proposed model under normal scenario.  [107] in B7-H3 knockdown cells grown in normoxic or perturbed conditions (hypoxia) for 24 h. (B) Higher lactate production has been shown by Li et al. [105] in Dppa4 overexpressed cancer cell, and (C) Similar altered behavior of lactate production has been depicted by our computational results in cancer cells compared to that in normal and perturbed ones.
We have further compared the model with a previous simulation result [71] for an environment in cancer cells, particularly, hypoxia condition [98]. Here, we have perturbed the present model to incorporate the mutation exhibiting hypoxia condition. Moreover, the simulation results have been compared with CE-MS measurements during hypoxia in human erythrocytes [71]. As hypoxia in human erythrocytes enhances the expression of certain glycolytic enzymes including hexokinase, aldolase and pyruvate kinase [71], we have applied these enzyme expression values in (0, 1) to the model following their activities as described in the previous investigation [71]. Fig. 6 depicts that our simulation results almost follow not only the CE-MS measurements but also the simulation results obtained by the study [71] in response to the enzyme activities due to hypoxia in human erythrocytes. We have calculated mean squared error (MSE) between our simulated concentrations of eight metabolites (G6P, fructose-6P, fructose 1, 6 bisphosphate, Dhap, 3PG, PEP, pyruvate, and lactate) with time, and the corresponding CE-MS measurements of the same due to aforementioned enzyme activities during hypoxia in human erythrocytes. The MSE values, corresponding to these eight metabolites, are (0.07337, 0.06031, 0.074389, 0.12327, 0.13879, 0.14829, 0.34182, and 0.05727). In addition, we have calculated the MSE between the simulation results of Kinoshita et al. [71]  In comparison with some other mathematical models in analyzing energy management in mammalian cancer cells, we have found that a constraint-based flux balance model [66] of CCM pathways in cancer cells has found a set of enzymes, such as lactate dehydrogenase, playing important roles in cancer growth. Similar results have been found using our model. We have elaborately discussed the results in the following sections. However, such a flux balance analysis based modeling does not depict the transient behavior of the molecules involved in the cancer pathway. Another differential equation based model [67] has explored a quantitative relationship between the hypoxia intensity and the intracellular lactate levels in cancer cells. Besides, it has predicted some important regulators of the glycolysis pathway only. Similar type of model [68] has been developed for targeting energy metabolism in pancreatic cancer.
Another model [69] based on optimization technique has shown the lactate-driven coupling for fulfilling energy requirements in cancer cells. However, unlike the proposed model, they do not consider three kinds of pathways with different response time. Moreover, they do not explore the effects of activating/deactivating the probable drug targets to control energy metabolism in cancer cells.
For more validation of the proposed model, we have compared our simulation results with in vitro observations from cultured human sertoli cells during insulin deprivation [72]. Such a condition in human sertoli cells reduces the expression level of lactate dehydrogenase (LDH). Besides, the expression level of glut1 has increased [72]. Incorporation of these two enzyme expressions in (0, 1) into our model, we can notice similar behavior of initial glucose consumption, pyruvate consumption and lactate production with the corresponding behavior in cultured human sertoli cells during insulin deprivation (Fig. 7).
In the context of present investigation, we have finally simulated the situation of knocking out certain enzymes, viz., pyruvate dehydrogenase, pyruvate carboxylase, acyl-CoA synthetase, fatty acid synthase, phosphoenolpyruvate carboxykinase 1 and succinyl CoA synthetase, to slow down oxidative phosphorylation through mitochondria leading to ''Warburg effect" [24,25]. In order to knock out these enzymes, we have considered the kinetic rate constants as zero for the reactions catalyzed by the above enzymes. Considering the aforesaid perturbed condition, we have monitored the behavioral alteration of different molecules compared to that in normal situation as depicted in Figures S3 and S4 in Supplementary material. Enhancement of pyruvate and lactate production as well as slowing down the production of FFA, citrate and succinate compared to that under normal condition indicates low oxidative phosphorylation. Subsequently, consumption of glucose increases to compensate less ATP production, which leads to reduction of glucose production compared to that under normal condition. In spite of this situation, ATP production decreases significantly. As most of the glucose is utilized through glycolytic flux, NADPH production decreases compared to that in normal mammalian cells. However, production of fructose 6P, acetyl CoA, ribose 5P, fructose 2,6 bisphosphate, PEP and protein kinase B (AKT) does not show any significant difference compared to that under normal condition. Interestingly, expression level of p53 increases to inhibit glucose transporter 1 (glut1) [99], which is responsible for transportation of glucose across the plasma membranes of mammalian cells [100]. Thus, higher expression of p53 leads to suppressing the production of glucose in glycolysis pathway. As a Table 4 Illustrating the significance of certain rational drug targets in terms of management of energy and cell proliferation in mammalian cancer cells.

Drug target Energy (ATP) production
Cell proliferation indicated by the production of ribose 5P and NADPH result, glucose concentration decreases rapidly (Fig. S3 in Supplementary material) leading to cell apoptosis [101]. Nevertheless, mammalian cancer cells somehow manage both energy, in the form of ATP, and cell proliferation to survive [102]. In this study, we aim at investigating the probable mutated regulation that drives mammalian cancer cells to survive despite ''Warburg effect". In this context, we have considered the aforesaid perturbed state space model of integrated biochemical pathways related to carbon metabolism as a MIMO plant of our interest. Thereafter, we have applied the GA controller to achieve high concentrations of ATP (0.95) and ribose 5P (0.70). Here, high ATP will supply constant energy in mammalian cancer cells to survive, while ribose 5P plays an important role in nucleotide synthesis promoting cell growth and proliferation [27]. Thus, the proposed MIMO plant along with GA controller mimics the altered nonlinear dynamics of mammalian cancer cells. We are now going to describe the probable mutated regulations of cancer cells in the following subsection.

Analysis of regulations in mammalian cancer cells
Here, we have considered the integrated CCM pathway with incorporation of possible mutation responsible for energy man-agement in cancer cells. The necessary pathway and information have been obtained from KEGG database [23] and literature. As a result, oncogenic somatic and germline mutation have automatically been taken care. The simulation results have found that ATP and ribose 5P production in mammalian cancer cells becomes sufficiently higher than in normal ones following the reference outputs to GA controller (Figs. 8A and 8D). In this situation, we have monitored the altered behavior of other molecules that assist cancer cells to survive and grow. We have found that the expression level of PHD decreases significantly (Fig. 9E) compared to that in normal cells. As PHD accelerates the degradation of hypoxia-inducible factor-1a (HIF-1a) [108,78], the expression level of HIF-1a is quite higher (Fig. 10B) due to less expressed PHD.
Subsequently, the expression levels of phosphatidylinositol-4,5-bisphosphate 3-kinase (PI3K) (Fig. 11D), AKT (Fig. 12B), mammalian target of rapamycin (mTOR) (Fig. 11F), MYC (Fig. 13C), and extracellular signal-regulated kinases (ERK) (Fig. 9F) increase significantly in cancer cells compared to that in normal ones. Besides, signal transducer and activator of transcription 3 (STAT3) and nuclear factor kappa-lightchain-enhancer of activated B cells (NF-jB) are highly expressed (Fig. 14F and 11B). On the other hand, p53 (Fig. 15B) shows under expression compared to normal mammalian cells. In support of our simulation results, we have found some previous investigations, including in vivo and in vitro experiments [28][29][30][31][32][33][34][35][36][37][38][39][40][41][42][43][44]64,65], which show similar kinds of behavior. Besides, we have collected the recent in vivo/in vitro data [54,56,57,55,[58][59][60][61][62][63] of human cancer cells to compare different mutations verified in experimental laboratories with our simulation results. Here, we have normalized the experimental data in ½0; 1. According to a clinical study [28], the data, collected from patients with breast, cervical and endometrial cancers at early stage, show high mortality rate of the patients having tumors with over expressed HIF-1a transcription factor. PHD inhibition promotes more stabilization of HIF-1a in cancer cells than in normal ones [29,30]. In this context, Kwon et al. [54] has found higher expression level of HIF-1a measured in human hepatocellular carcinoma (HCC) specimens compared to the non-cancerous tissue specimens (Fig. 10A). Previous studies [31][32][33] have shown that PI3K, AKT and mTOR express significantly higher in tumor tissues of patients suffering from ovarian, gastric and prostate cancer compared to normal individuals. Experiment [55] with human colon cancer sample has depicted higher expression of PI3K compared to normal ones (Fig. 11C). Immunohistochemical scores measured by Mao et al. [56] has depicted that AKT expression is significantly higher in pancreatic cancer (91 cases) than that in normal pancreas (51 cases) (Fig. 12A). Besides, Roy et al. [57] has shown the deferentially expressed AKT isoforms in normal and malignant oral tissues through immunohistochemical analysis of the human samples (Fig. 12A). On the other hand, higher mTOR expression has been found in colon cancer compared to normal one as per the in vivo/in vitro experiment performed by Zhang et al. [55] (Fig. 11E).
Evidences [34,35,65] suggest that over expression of MYC triggers certain genes to promote growth and proliferation of cancer cells. Similar behavior of MYC in cancer cells has been found in two recent in vivo/in vitro experiments [58,59] (Fig. 13A and  13B). It is reported that tumor suppressor p53 is under expressed [36] in cancer cells due to lysine methylation [37]. Moreover, evidences [38,39,40] have shown that highly expressed ERK in cancer cells promotes over expression of N-cadher protein involved in metastasis. Subsequently, previous investigations [41,42] have demonstrated that enhanced expression of STAT3 signaling protein either inhibits apoptosis or accelerates cell proliferation, angiogenesis, and finally metastasis. Consequently, it leads to cancer initiation and progression. Besides, over expression of NF-jB plays an important role in this context [39,43,44]. The previous in vivo/ in vitro studies have supported the aforementioned claim by illustrating over expressed STAT3 in certain cancer cells, such as breast cancer [60], gastric cancer [61] and colon cancer [62] cells compared to normal ones (Fig. 14E). Besides, the evidence [63] has also depicted that relative expression of NF-jB significantly increases in cancer cells (CT-2A astrocytoma) compared to normal ones (Fig. 11A). The present simulation results demonstrate higher expression of hexokinase (HK) and glut1 (Figs. 16D and F) in cancer cells. Previous in vivo/in vitro result [106] depicting higher expression of HK mRNA levels (Fig. 16C) in cancer compared to normal liver tissue has validated the computational result. This experiment [106] has also shown higher glut1 mRNA levels (Fig. 16E) as 92-fold higher in Meta specimens compared to normal liver tissue. Here, enhanced expressions of HIF-1a, PI3K, AKT, MYC and mTOR are responsible for activating HK and glut1 more in cancer cells. As a result, glucose consumption increases (Fig. 17C). Consequently, glucose production decreases as depicted in Fig. 15E. Previous investigations [27,45] support our claim. They have claimed that enhanced glut1 increases the utilization of glucose by anabolic pathways. Besides, highly expressed MYC activates lactate dehydrogenase (LDH) more in cancer cells. It leads to fermentation of glucose [109,25] through enhanced glycolysis [46]. Consequently, lactate production increases. The present results have followed the above claims by showing higher expression of LDH (Fig. 16B) as well as enhanced lactate production (Fig. 18C). In this context, Lim et al. [107] has measured higher glucose uptake (Fig. 17A) in B7-H3 knockdown cells grown in normoxic or perturbed conditions (hypoxia) for 24 h. Besides, Li et al. [105] has shown that glucose uptake increases significantly in Dppa4 overexpressed cancer cells (Fig. 17B). This experiment [105] has also shown that LDH has significantly been up regulated (Fig. 16A) in Dppa4 overexpressed cancer cells. As a result, higher lactate production (Fig. 18B) has been found here [105]. Lim et al. [107] has also found similar higher lactate production (Fig. 18A) in B7-H3 knockdown cells as mentioned before. These experimental results provide a strong support to the present computational results.
Higher expression levels of phosphofructokinase 1 (PFK1) (Fig. 14B), phosphofructokinase 2 (PFK2) (Fig. 9A), and glyceraldehyde-3-phosphate dehydrogenase (Fig. 15A) also confirm enhanced glycolysis in mammalian cancer cells. In support of the computational results, a previous in vivo/in vitro result [104] has depicted higher PFK1 expression (Fig. 14A) in breast cancer and paracancer tissues, expressed as units per gram of protein (U/gprot). Consequently, higher amount of fructose 6P is utilized to produce more fructose 2,6 bisphosphate and fructose 1,6 bisphosphate than in normal cells. Besides, higher amount of fructose 2,6 bisphosphate accelerates the break down of fructose 6P into fructose 1,6 bisphosphate [96]. That is why simulation result has shown reduction of fructose 6P production (Fig. 8B). Similarly, the concentrations of PEP drops (Fig. 8E) due to its over consumption to produce higher amount of pyruvate and ATP (Figs. 8F and 8A) compared to normal cells. Here, expression of pyruvate kinase switches alternatively from low to high and vice versa as depicted in Fig. 14D. Evidences [47][48][49]64] have shown that pyruvate kinase (M2 isoform) switches to its inactive dimer form or active tetrameric form according to the requirements of mammalian cancer cells. When pyruvate kinase (M2 isoform) is in tetrameric form, the flux through glycolysis enhances with sufficient amount of ATP production. As a result, production of intermediate glycolytic metabolites, such as glyceraldehyde 3P (Fig. 8C) and fructose 1,6 bisphosphate, increases. Conversely, dimer form of pyruvate kinase (M2 isoform) promotes macromolecular synthesis from glycolytic intermediate metabolites to continue cell growth and proliferation. An evidence [104] has also reported about higher pyruvate kinase expression (Fig. 14C) in breast cancer and paracancer tissues, expressed as units per gram of protein (U/gprot). Thus, cancer cells manage energy in the form of ATP to survive in spite of slow oxidative phosphorylation under consideration. Here, reduction of reduced nicotinamide adenine dinucleotide (NADH) production ( Fig. 10D) indicates successful incorporation of slow oxidative phosphorylation [50] into the proposed model. An in vivo/in vitro experiment, performed by Sumi et al. [103] depicting lower concentration of NADH (Fig. 10C) in cancer cells compared to normar ones, supports the computational result. According to the simulation results, the enzymes (proteins) catalyzing PPP, such as glucose-6-phosphate dehydrogenase (Fig. 9D), phospho-gluco dehydrogenase (Fig. 15C), ribose 5P isomerase (Fig. 9C), transketolase and transaldolase ( Fig. 15D and 9B), have shown over expression in cancer cells for macromolecular precursors required for cell growth and proliferation. Under expressed p53 plays an important role in this context. High concentration of ribose 5P (Fig. 8D) confirms enhanced production of cell building material, including DNA, RNA, nucleic acids and histidine, in mammalian cancer cells. Evidences [51][52][53] have shown that relatively higher expressions of glucose-6-phosphate dehydrogenase, phospho-gluco dehydrogenase, ribose 5P isomerase, transketolase and transaldolase help cancer cells in generating high amount of NADPH and ribose 5P, which are responsible for reactive oxygen species (ROS) reduction as well as the production of high levels of nucleotides for DNA synthesis and repair. It leads to resistance against certain cancer therapies resulting in enhancement of oxidative stress or DNA damage. In this context, it should be mentioned that higher amount of ROS can frequently be observed in cancer cells helping in activation of oncogenes and metastasis. However, further enhancement of ROS beyond a certain threshold induces cell death [110]. Finally, we have summarized the altered regulations of different metabolites, transcription factors and genes in cancer cells compared to that in normal ones in Table 3.

Analysis of certain rational drug targets in terms of management of energy and cell proliferation in mammalian cancer cells
Here, we are going to discuss the effects of certain rational drug targets on cancer cells in terms of management of energy and cell proliferation from simulation point of view. These results may have significant implications during in vivo and/or in vitro experiments. When the proposed model computationally meets the reference target concentrations of ATP and ribose 5P for GA controller, we have observed that the model has mimicked the altered regulation of cancer cells as previously discussed. At that moment, we have set reference target expression levels (high or low) of certain proteins/enzymes as drug targets using GA controller. Although the model is quite capable of analyzing the effect of any drug target, we have considered only six drug targets among others for our study to restrict the size of present article. Here, we have monitored the effects of deactivating pyruvate kinase (reference expression level 0.02), glucose-6-phosphate dehydrogenase (reference expression level 0.02), transketolase (reference expres- sion level 0.09), ribose 5P isomerase (reference expression level 0.09) and glucose-6-phosphate isomerase (reference expression level 0.03). Finally, the effect of activating pyruvate kinase (reference expression level 0.98) has also been observed and analyzed. Table 4 summarizes the significance of the drug targets under consideration in energy supply and cell proliferation.
Case 1 (Pyruvate kinase deactivation): According to the simulation results as depicted in Fig. 19, we have found that deactivation of pyruvate kinase cannot reduce the concentration of ribose 5P (Fig. 19A) and NADPH (Fig. 19C). In other words, high proliferation of mammalian cancer cells continues. Besides, enhanced NADPH still prevents from ROS production helping cancer cells to survive in spite of oxidative stress. Although the concentrations of ATP (Fig. 19A) and lactate (Fig. 19C) reduce at the beginning, after a while they increase again significantly. Even pyruvate kinase deactivation cannot reduce high glucose utilization (i.e., high consumption and less production) (Fig. 19B) in mammalian cancer cells. Moreover, high expression of glut1 and under expression of tumor suppression protein p53 (Fig. 19D) cannot be reversed in this case. Thus, pyruvate kinase deactivation may not be a good choice as a drug target for mammalian cancer therapy.
Case 2 (Glucose-6-phosphate dehydrogenase deactivation): Glucose-6-phosphate dehydrogenase is the key enzyme (protein) to utilize G6P through PPP. It leads to producing high amount of NADPH and ribose 5P so that nucleotides and fatty acid can sufficiently be synthesized to maintain intracellular redox homeostasis [112]. The present results (Fig. 20) have shown that the concentrations of ribose 5P (Fig. 20A) and NADPH ( Fig. 20C and 21B) decrease due to deactivation of glucose-6-phosphate dehydrogenase. In support of this result, a recent in vivo/in vitro study has shown how the relative activity (expression level) of glucose-6-phosphate dehydrogenase enzyme controls the NADPH production [111]. Lower expression of glucose-6-phosphate dehydrogenase leads to reduction of the concentration of NADPH (Fig. 21A). Here, expression of glucose-6-phosphate dehydrogenase and NADPH concentration have been assessed in A549/DDP cells according to various concentrations of 6-Aminonicotinamide (6-AN). Besides, expression level of p53 increases (Fig. 20D), whereas glut1 expression decreases (Fig. 20D). These results indicate reduction of cell growth and proliferation due to glucose-6phosphate dehydrogenase deactivation. However, high amount of glucose (Fig. 20B) is utilized through glycolysis to generate sufficient ATP (Fig. 20A) for cell survival. Subsequently, cell fer- mentation continues through high amount of lactate production (Fig. 20C). A previous study [113] has shown reduction of cell proliferation by silencing glucose-6-phosphate dehydrogenase in the human breast cancer cell line MCF7. Thus, deactivation of glucose-6-phosphate dehydrogenase as a possible drug target may be efficient to reduce cell growth and proliferation but not be effective in reduction of energy supply and fermentation. Case 3 (Transketolase deactivation): Transketolase has similar importance as glucose 6-phosphate dehydrogenase to maintain cell growth and proliferation leading to metastasis [114]. Simulation results (Fig. 22) demonstrate that the concentration of ribose 5P (Fig. 22A) decreases significantly in accordance with the deactivation of transketolase. Subsequently, enhanced expression of p53 (Fig. 22D) conveys that cell growth and proliferation are inhibited in this case. Reduction of NADPH concentration (Fig. 22C) signifies enhancement of ROS preventing from survival of cancer cells by oxidative stress. Moreover, ATP (Fig. 22A) and lactate (Fig. 22C) production also decrease through reduced glucose utilization (Fig. 22B). Besides, expression level of glut1 decreases as depicted in Fig. 22D. Previous investigations [113,115,116] have shown similar effects due to transketolase silencing. Thus, it is clear that targeting transketolase may be a potential choice for future cancer therapy. Case 4 (Ribose 5P isomerase deactivation): We have already discussed that up regulation of ribose 5P isomerase plays an important role in cell growth and proliferation of a cancer patient. A previous study [117] has claimed that expression level of ribose 5P isomerase is enhanced in colorectal cancer. Our simulation results (Fig. 23) have shown the effect of silencing ribose 5P isomerase in cancer cells. We have observed that enhanced p53 expression and decreased glut1 expression (Fig. 23D) indicate slowing down of proliferation and growth of cancer cells. In this context, reduction in NADPH (Fig. 23C) and ribose 5P (Fig. 23A) concentration confirm our claim. Subsequently, glucose utilization (Fig. 23B) and cell fermentation (i.e., lactate production as depicted in Fig. 23C) decrease. However, energy supply in the form of ATP (Fig. 23A) in cancer cells is somehow managed in spite of the deactivation of ribose 5P isomerase. Consequently, ribose 5P isomerase might be considered as a biomarker for targeted cancer therapy and prediction. Case 5 (Glucose-6-phosphate isomerase deactivation): Glucose-6-phosphate isomerase or phosphoglucose isomerase (PGI) is a glycolytic enzyme that directs G6P flux into the glycolysis branch. As a result, G6P breaks down into fructose 6P. Deactivation of Glucose-6-phosphate isomerase results in slowing down of glucose utilization (Fig. 24B) as well as lactate (Fig. 24C) and ATP (Fig. 24A) production. Thus ''Warburg effect" is somehow reversed. Besides, reduced glut1 expression level along with enhanced p53 expression (Fig. 24D) may slow down growth and proliferation rate of cancer cells. In this context, reduction of NADPH (Fig. 24C) and ribose 5P (Fig. 24A) concentration supports our claim. Here, a previous investigation [118] has shown that down regulation of Glucose-6-phosphate isomerase may suppress ''Warburg effect". Besides, oxidative phosphorylation is activated. However, its impact on tumor growth is minimal except in the case of hypoxia. Thus, glucose-6-phosphate isomerase can be a potential drug target for future cancer therapy. Case 6 (Pyruvate kinase activation): Activation of pyruvate kinase (Fig. 25) may slow down cell growth and proliferation as well as cell fermentation. Decreased concentrations of NADPH, lactate (Fig. 25C) and ribose 5P (Fig. 25A) confirm our claim. Besides, decreased glucose utilization (Fig. 25B) and glut1 expression (Fig. 25D) support the fact. Here, expression level of p53 (Fig. 25D) also increases to slow down cell growth and proliferation. However, ATP is produced sufficiently in this case as depicted in Fig. 25A. Evidence [47] has claimed that cancer cell growth and proliferation are inhibited by activation of pyruvate kinase (M2 isoform). Thus, pyruvate kinase (M2 isoform) activators may be considered as possible significant drugs for oncogenic treatment.

Conclusion
In this study, we have successfully integrated three types of biochemical pathways, viz., metabolic, signaling and gene regulatory networks, keeping their three timescale nature intact. Here, we have developed the integrated state equations considering appropriate timescales as well as all possible perturbations present in the contemplated integrated biochemical pathway. Besides, depending on the training dataset generated by solving the pathway ODEs, SVR based MIMO model has been developed. The MIMO model can mimic the transient nonlinear dynamic behavior of the integrated biochemical pathway under consideration. Moreover, with the help of the GA controller, the model can predict the effect of drug targets applied to complex diseased cells. In order to investigate the effectiveness of the model, we have used our model to explore how mammalian cancer cells are able to manage their growth, proliferation and energy supply to survive. In this context, ''Warburg effect" [24,25] has been taken into account. The simulation results have depicted that the model has not only captured the key regulations, but also has been able to predict certain possible drug effects in terms of energy and cell proliferation management in mammalian cancer cells.
According to the results, the proteins or genes HIF-1a, HK, glut1, AKT, glyceraldehyde-3-phosphate dehydrogenase, phospho-gluco dehydrogenase, ERK, ribose 5P isomerase, mTOR, glucose-6phosphate dehydrogenase, STAT3, NF-jB, PI3K, MYC, LDH, PFK1, PFK2, transketolase and transaldolase are up regulated in cancer cells. Besides, PHD and p53 are down regulated. Switching of pyruvate kinase (M2 isoform) between its two oligomeric form, viz., inactive dimer and active tetramer, plays an important role in managing proliferation, growth and energy in mammalian cancer cells. These results have been validated through previous investigations involving in vivo and in vitro experiments [24,27,26,28-4 5,25,46-54,56,57,55,58-65]. Besides, other mathematical models [66][67][68][69] support the results derived by the proposed model. Among six drug targets under consideration, deactivation of transketolase and glucose-6-phosphate isomerase may be the most potential to slow down cancer progression by reducing cell proliferation, growth, fermentation and energy supply. On the other hand, pyruvate kinase (M2 isoform) activation and ribose 5P isomerase deactivation may reduce cell growth, proliferation and fermentation during cancer. However, they may not be able to stop energy supply in mammalian cancer cells. Although deactivation of glucose-6phosphate dehydrogenase may slow down cell growth and proliferation, it may fail to stop fermentation and energy supply in the malignant cells. In this context, pyruvate kinase deactivation may be an awful choice as a rational drug target for future cancer therapy.
Finally, it may be mentioned that the experimental values of different kinetic parameters are still unknown or poorly documented. In this study, we have estimated these values by trial and error based on previous knowledgebase. Some of the parameter values considered here have been seen to be close enough with those estimated by the method of Lillacci et al. [70]. The other parameter values could not be checked due to unavailability of required experimental observations. However, trial and error based technique is time consuming and difficult to perform because of large number of such kinetic parameters. This is a drawback of the present methodology. Thus more in vivo/in vitro parameter values are needed, and/or appropriate parameter estimation methods based on the theory of machine learning and control theory can be developed for designing more accurate MIMO plant for an integrated biochemical pathway.