Integrative Computational Framework for Understanding Metabolic Modulation in Leishmania

Background: The integration of computational and mathematical approaches is used to provide a key insight into the biological systems. Through systems biology approaches we seek to find detailed and more robust information on Leishmanial metabolic network. Forman/Forman-Ricci curvature measures were applied to identify important nodes in the network(s). This was followed by flux balance analysis (FBA) to decipher important drug targets. Results: Our results revealed several key high curvature nodes (metabolites) belonging to common yet crucial metabolic networks, thus, maintaining the integrity of the network which signifies its robustness. Further analysis revealed the presence of some of these metabolites, MGO, in redox metabolism of the parasite. Being a component in the glyoxalase pathway and highly cytotoxic, we further attempted to study the outcome of the deletion of the key enzyme (GLOI) mainly involved in the neutralization of MGO by utilizing FBA. The model and the objective function kept as simple as possible demonstrated an interesting emergent behavior. The non-functional GLOI in the model contributed to “zero” flux which signifies the key role of GLOI as a rate limiting enzyme. This has led to several fold increase production of MGO, thereby, causing an increased level of MGO•− generation. Conclusions: The integrated computational approaches have deciphered GLOI as a potential target both from curvature measures as well as FBA which could further be explored for kinetic modeling by implying various redox-dependent constraints on the model. Furthermore, a constraint-based FBA on a larger model could further be explored to get broader picture to understand the exact underlying mechanisms. Designing various in vitro experimental perspectives could churn the therapeutic importance of GLOI.


INTRODUCTION
Metabolic networks are becoming of much interest amongst the researchers due to their unique feature of connecting every node (the metabolites) with the links (the reactions) that are catalyzed by specific gene products responsible for growth and maintenance of a cell (Mahadevan and Palson, 2005). Metabolic networks can be used to shed light on finding various disease mechanisms via the identification of essential genes by implying perturbations on networks, thereby, characterizing the structure-function relationships (Watts and Strogatz, 1998;Barabási and Oltvai, 2004;Dorogovtsev and Mendes, 2013). Structure-function relationship of a network in terms of network robustness and reliability is strongly linked to its geometric and topological properties and can be studied by using graph theory (Barabási and Oltvai, 2004). Recent technological advances in graph theory have allowed us to analyze and describe highly complex systems at smaller and more detailed scales and have been utilized to study global impact of long term sickle cell disease on brain (Case et al., 2018), to diagnose pre-symptomatic Alzheimer's disease (Bytautiene, 2013), to understand the human brain network (Vecchio et al., 2017). Recently, the development of geometry based measures is significantly attracting the focus of researchers to characterize the structural aspect of complex networks.

Curvature-Based Methods for Complex Networks
In geometry based measures of complex networks, curvature is one of the notions that are being explored to understand the complexity of graphs. Curvature plays central role in Riemannian geometry since it represents a measure to quantify the deviation of a geometrical object from being flat (Tannenbaum et al., 2016). Curvature measures can be used to quantify the robustness and thereby the functionality of networks.
Among several types of curvature notions, Ricci curvature is known to be the most useful for analyzing the complex networks (Burago et al., 2001;Saucan and Appleboim, 2005;Ollivier, 2009Ollivier, , 2010Ollivier, , 2011Ni et al., 2015;Sandhu et al., 2015;Sreejith et al., 2016). Ricci curvature measures deviance of geodesics (shortest path) relative to Euclidean shortest-paths and is related to mass transport or entropy (Wasserstein metrics) (Evans, 2001;Villani, 2003). High Ricci curvature is typically found near network hubs, where an anchor node (high curvature node) is connected to many existing nodes in close proximity. Principally, curvature is very close to robustness; therefore, removal of high curvature node/edge will result in collapse of network.
Ollivier-Ricci (Ollivier, 2009(Ollivier, , 2011 and Forman-Ricci (Forman, 2003), two different discretization of Ricci curvature, give efficient solutions for network geometrizations. In case of undirected networks, Ollivier-Ricci curvature has proved its applicability in various network analyses (Loisel and Romon, 2014;Ni et al., 2015;Sandhu et al., 2015;Gao et al., 2016). On the other hand, Forman-Ricci curvature has been introduced as a tool for undirected (Sreejith et al., 2016(Sreejith et al., , 2017a as well as directed (Weber et al., 2017a,b) network analyses. Irrespective of their applicability and handling the large networks, the two curvature notions demonstrated high correlation (Samal et al., 2018). However, Forman-Ricci curvature is a faster computation method and can be utilized in larger real-networks. Current study exploits the application of Forman-Ricci curvature and systems biology to capture the behavior of Leishmanial metabolic network.

Leishmaniasis and Systems-Biology
Neglected infectious diseases have affected at least a billion of human populations worldwide (Narain et al., 2010;WHO, 2011) and are of primary concern due to the lack of effective and affordable drug regimens. Among them, cutaneous leishmaniasis (CL), a very common clinical form of leishmaniasis, has always been overlooked as a major public health problem due to its non-fatality. The causative agent of CL, a protozoan parasite, Leishmania major has a digenetic lifecycle and lives in two hosts, sandfly, Phlebotomus argentipes, and human, in the form of flagellated promastigotes (procyclic phase) and non-flagellated amastigotes (metacyclic phase), respectively.
When inside the human host, the parasite has to undergo the oxidative stress generated by host macrophages (Figure 1). This stress is also attained from free radical generation within the parasite in many ways mainly via electron transport chain. Moreover, there are several highly reactive metabolites such as glyoxals that contribute to the generation of free radicals. These metabolites can react with the other metabolites and enzymes and not only cause inactivation of enzymes but also the formation of free radicals. To deal with these types of internal and external trauma the parasite has evolved with its unique redox machinery that helps in neutralization of free radicals. The main component of this antioxidant defense mechanism are trypanothione synthase (TryS), trypanothione reductase (TryR), tryparedoxin peroxidase (TryP), tryparedoxindependent peroxidase (TDPx), and tryparedoxin (TXN). The center of these components is a central reductant, T[SH] 2 that regulates the activation of these enzymes (Krauth-Siegel and Inhoff, 2003).
Due to the incapability and increasing drug resistance of the current therapy regime techniques, the necessity to develop novel, more efficient and affordable anti-leishmanial drugs along with potential molecular drug targets with therapeutic value are needed. The difference between Leishmania and human host redox metabolism attracts researcher's interest to explore this mechanism for discovering novel drug targets and design potential inhibitors against the redox enzymes without crossing the host machinery (Ascenzi et al., 2003).
To study the role and effect of other intermediates in the redox network, systems biology approaches are nowadays being used (Kabra et al., 2018). Systems biology approaches have already been utilized to study systems dynamics of complex redox metabolic pathways in many organisms (Hädicke et al., 2011;Pillay et al., 2013;Wang et al., 2018). Before implementing many available network analyses techniques like, Genome scale metabolic modeling (GEMs) (Hädicke et al., 2011;Pillay et al., 2013), constraint based methods, kinetic pathway modeling (Ho, 2008;Resat et al., 2009;Kumar et al., 2019), the knowledge and availability of all possible components and their interconnections is a pre-requisite. Recently, kinetic modeling and constraint-based metabolic modeling have gained researcher's focus. In case of kinetic modeling, the availability of kinetic data and rate laws is required. However, due to the lack of sufficient kinetic data the use of metabolic kinetic modeling on larger scale networks is limited. On the other hand, constraint-based analysis such as flux-based analysis (FBA) requires only the knowledge of stoichiometry of the metabolites and can be easily applied to even larger scale metabolic networks. FBA calculations rely on the assumption of steady-state growth and mass balance and can be used to predict the flux of a metabolite flowing through a metabolic network. To predict the growth rate or to find a specific cellular function, FBA tries to find possible solutions to optimize the stated objective function(s) and gives quantitative insights into the genotype-phenotype relationship of the model. Although, FBA is incapable of predicting the concentrations of metabolites and does not account for any regulatory effects like activation/inactivation of enzymes or gene-expression regulations, there are extensive embodiment of literature showing its applications in many fields including drug target identification (Chavali et al., 2008;Fatumo et al., 2009;O'Brien et al., 2015;Sharma et al., 2017;Subramanian and Sarkar, 2017;Tewari et al., 2017).
In current work, we attempted to explore the application of Forman and Forman-Ricci curvature on metabolic network of Leishmania. Using curvature measures, we tried to identify important metabolic pathways contributing to the overall network robustness by recognizing "important" metabolic hubs in the network. Through this, we seek to spot metabolites and enzymes having essential role in the survival of the parasite. Further, to predict the effect of the absence of these "important" metabolites, constraint-based FBA is used to uncover the importance of enzymes and their probable role as drug target for therapeutic significance (Figure 2). Also, our systems pharmacology modeling on the kinetic model of the selected drug target using known inhibitors may provide significant insight for further clinical interventions.

Data Sets and Model Reconstruction
Initially, Leishmania major metabolic network iAC560 (Chavali et al., 2008) was retrieved from BioModels Database and imported to Cell Designer (4.4) (Funahashi et al., 2003). The model consisted of 1112 reactions and 1147 metabolites. This model was used to construct a bipartite directed graph that considered as master model. The master model was further used for constructing small scale models (Figure 2).

Generation of Connectivity Matrix
Connectivity matrices for master and other models were generated by considering the connections between metabolites (substrate to product). A value of "1" was assigned for connection, if a metabolite (substrate) is connected (producing or converting into) to other metabolite (product), and "0" for no connection.

Network Topology
The reconstructed bipartite directed graphs were visualized in Cytoscape (3.6.0) (Shannon et al., 2003) and their topological information (network properties) of nodes and edges were obtained using graph theory followed by calculation of various network properties. We used Forman curvature and Forman-Ricci curvature for the analysis of our bipartite directed graphs. For any given directed bipartite graph we calculated (Samal et al., 2018): Where, F is the Forman curvature of the directed edge e = − − → v 1 v 2 that originates from node v 1 and terminates at node v 2 . Only those directed edges were considered for calculation that either terminate at node v 1 or originate at node v 2 . w e denotes the weight of the edge e under consideration. w v1 and w v2 denote the weights associated with the nodes v 1 and v 2 , respectively. e v1 ∼ e and e v2 ∼ e denote the set of edges incident on nodes v 1 and v 2 , respectively, after excluding the edge e under consideration which connects the two nodes v 1 and v 2 . Furthermore, self-loops or self-edges on nodes v 1 and v 2 , edges facing opposite direction, were completely ignored.

Forman Curvature on Networks
To distinguish between incoming, E I,v , and outgoing, E O,v , edges for a given node v, unnormalized Forman curvature (the sum of the curvature of all edges incident or outgoing on that node), In Forman curvature F I(v) and Out Forman curvature F O(v) were computed as follows: Frontiers in Bioengineering and Biotechnology | www.frontiersin.org And, the total flow on a given node was obtained as:

Forman-Ricci Curvature on Networks
Moreover, for determining Forman-Ricci curvature, FR(v) (normalized Forman curvature) (Sreejith et al., 2017a,b) of a node v, the sum of the curvature of all edges incident or outgoing on that node was divided by the degree of that node: Where, FR I(v) and FR O(v) are In and out Forman-Ricci curvature, respectively, deg(v) is in and out-degree and F(e v ) is the Forman curvature of edge e v , and e v ∼ v represents the set of edges incoming or outgoing on that node v. Subsequently, total flow on node v was obtained as follows: Flux Balance Analysis

Reconstruction of the stoichiometric model
For the reconstruction of stoichiometric model COBRA Toolbox (Vlassis et al., 2014;Heirendt et al., 2018) was used. Necessary exchange reactions and transport reactions for uptake and excretion of typical substrates and products, respectively, were introduced. The network was devoid of any organism specific internal compartments and only external and internal environment were taken into account. Stoichiometric model is represented in the form of stoichiometric matrix, S, that is a mathematical representation of the network. In this matrix, each reaction is a column and each metabolite is a row. And the participation of the metabolites in the corresponding reactions is denoted by their stoichiometric coefficients where "−1" is used for the consumption and "+1" for the production of the metabolite. In our model, the upper and lower bounds of reversible reactions were set between +1,000 and −1,000, and for irreversible reaction between 0 and +1,000. Uptake and excretion reactions were also set between −1,000 and +1,000. The upper and lower bounds for any reaction define the total space allowed for flux distribution (Orth et al., 2010).

FBA simulation
Once the model is reconstructed, constraint-based FBA simulation was performed in COBRA Toolbox in MATLAB (Becker et al., 2007). FBA assumes steady-state kinetics and uses linear programming (LP) based optimization to determine the flux distribution to solve a given metabolic objective function by maximizing or minimizing it (Orth et al., 2010). Hence, the LP problem was formulated as, Maximize Z, Subject to Where, Z is the objective function to be maximized (or minimized), S is the stoichiometric matrix of m x n, v represents the flux vector that is controlled through enzyme capacity constraints v min and v max representing lower and upper bounds, respectively. After conversion into mathematical form, the simulation was performed to maximize the objective function. Initially the objective function was carrying zero flux indicating metabolic gaps in the model. The required demand reactions were added to the network to refine the model till the non-zero value of the objective function is attained.

Systems Pharmacology Modeling
The reactions for the model were input in COPASI v4.19.140 (Hoops et al., 2006) followed by assignment of appropriate kinetic laws. All the reactions were considered as irreversible and the reversible reactions were broken into two reactions for maintaining uniformity in the kinetic model ( Table S1). The concentrations, time units and reaction fluxes in the model are in molar (M), seconds and M/s, respectively. All the kinetic parameters were obtained and/or calculated from published literature (Data Sheets 1, 2, Table S2). To observe the effect of different parameters on the models sensitivity analyses was performed using deterministic LSODA ODE solver (Hindmarsh, 1983).

Model Reconstruction and Network Topology
A directed master network, M, was built after editing iAC560 by doing following: all transport and exchange reactions were removed; all self-loops, currency metabolites (ATP, ADP, CO 2 , etc.), electron carriers (NADP + , NADPH, NAD, NADH), Waters (H 2 O), and electron transporters (H + , PO +3 4 , HCO − 3 ) were removed; each reversible reaction in the network was converted to two irreversible backward and forward reactions; then, each irreversible reaction was considered as directed edge connecting the substrate metabolite to the reaction or reaction to the product metabolite; duplicate reactions and metabolites were carefully observed in the network and removed. Moreover, each directed edge was assigned with a weight corresponding to the stoichiometry of the involved metabolite (substrate or product) in the reaction under consideration. This model consisted of 226 metabolites (Table 1, Figure 3, Table S5).
Further, to get a deeper insight of the Leishmanial metabolic network, the master model, M, was subdivided into small directed bipartite graphs (Table 1, Figure S1). The network topology of networks is summarized in Table 1 and Table S3. The clustering coefficient (ranges between 0 to 1; no to high interconnectivity) of the master network, M, was 0.038. Among the models M1-M7, only M7 had the lowest clustering coefficient 0.01 ( Table 1). The clustering coefficient of M7 indicates that the network has lesser connectivity in comparison to other networks, showing the transmission of the signal is linear and in one direction.
Further, betweenness centrality calculated for each node in M network pointed out important metabolites. Nodes having high betweenness centrality are known to act as bridge between nodes for transmission of information in the network. In our M model several of these metabolites (HTA, MG, Pyr, T[SH] 2 ) with high betweenness centrality were part of M7 model (Table S3) suggesting the importance of M7 model .

Curvature Distribution in Master Network
Normalized (Forman) and unnormalized (Forman-Ricci) Forman curvatures of directed edges and nodes were calculated in all reconstructed models (Table S4). Forman curvature of nodes and edges in network "M" was computed and their distribution was plotted to analyze their frequency (Figure 4). From the plots, it was observed that most of the nodes and edges   were comprised of negative curvature values. Although, the distribution of Forman curvature of both nodes and edges was broad, the peaks with high frequency were concentrated toward (or at) "zero." The distribution of Fi and Fo of nodes were found to be slightly different than Forman curvature of nodes however the peaks between 0 and −2 were found to be with frequency >50. Similar was the case with Fi and Fo of edges where one peak in both the plots had frequency >50. Further extraction of these peaks demonstrated that metabolites belonging to amino acid metabolism, sterol synthesis, fatty acid biosynthesis and degradation, glycerolipid and glycerophospholipid metabolism pathways were part of these high peaks. Interestingly, peaks appearing far from "zero, " those with higher Forman curvature, consisted of edges and nodes related to glycolysis, MG metabolism and T[SH] 2 metabolism pathways. Furthermore, these metabolites were also spotted as important vertices from Forman curvature plot of nodes in "M" model ( Figure 5)   (Table S4).

Correlation Between Forman Curvature and Common Network Measures
Various common statistical topological measures were calculated for all networks and studied to find any significant correlation with the curvature measures ( Table 2).

Master Model
The correlation between Forman curvatures and betweenness centrality, closeness centrality, In-degree and Out-degree were computed in "M" model (Table 2, Figure 6). A high negative correlation was observed for Fi and FRi with In-degree  ( Figures 6B,F). Moderate negative correlation was reported for Fi, FRi, Fo with betweenness centrality (Figures 6A,C,E). Interestingly, In-degree and Out-degree have also shown moderate correlation with betweenness centrality but the R-value was shifted to positive side (Table 2). Moreover, in contrast, Fo and FRo have shown very less or near negligible correlation with Out-degree and betweenness centrality, respectively (Figures 6D,G,H). Forman curvature of nodes was also negatively correlated with high magnitude to betweenness centrality than that of closeness centrality suggesting the importance of the former in finding important nodes in the network ( Table 2, Table S4).

Small Scale Models
It was noted that, in all models (M1-M7) except M3 (−0.576), the correlation between Forman curvature of nodes and clustering coefficient was very weak negative in comparison to correlation between Forman curvature and degree. Interestingly, the highest negative correlation between Forman curvature and degree was mainly noted for model M5, M6 and M7 (−0.854, −0.782, and −0.903, respectively) ( Table 2). Further, weak negative correlation coefficient between closeness centrality and Forman curvature of nodes was obtained.
On the other hand, a good negative correlation of betweenness centrality with Fi and FRi was observed only in case of M5, M7 and M4, M7, respectively. Surprisingly, in contrast to Fi and FRi, correlation for Fo and FRo with out-degree and betweenness centrality was noted very weak to moderate negative. Among the small scale models, only M7 indicated high negative correlation for Fo and FRo with betweenness centrality. It was found similar to the correlation for betweenness centrality and In-and Outdegree for M7 model (high positive correlation) ( Table S4).

Robustness of Models
Forman and Forman-Ricci curvatures across all edges and nodes of the network were computed. High curvature nodes are said to be the backbone of the network those act as bridges between major network communities (Varma and Palsson, 1993) (Figures 6,  7, Figure S1). Hence, for investigating the robustness in terms of "network connectivity, " the nodes were removed on the basis of decreasing strength in terms of Forman curvature, Forman-Ricci curvature and clustering coefficient. Nodes with low curvature and low clustering coefficient were removed first. Highly connected network is represented as "1, " and a higher value is an indication of disconnected network. In Figure 8, we plotted network connectivity against the nodes removed on the basis of Forman curvature and clustering coefficient. It was observed that removal of high curvature nodes caused faster disintegration of the network. Nodes deletion according to Forman-Ricci curvature gave similar result as with Forman curvature. However, nodes removal on the basis of clustering coefficient has also caused network destruction but slower than Forman curvature (Figure 8) (Sreejith et al., 2016).
In our master model we found more than 25 nodes with high curvature values (Table S3). We noted that most of the nodes from these were related to "M7" model that consists of glycolysis, T[SH] 2 synthesis, and methylglyoxal metabolism pathways, signifying the importance of "M7" model .

Characteristics of the Reconstructed Network
As stated in previous section, the role of several metabolites including MGO and T[SH] 2 (part of M7 model) via Forman curvature and Forman-Ricci curvature measures was found important for maintaining the integrity of the master model (Table S3). Since our focus for performing FBA was to investigate the importance of glyoxalase pathway in redox metabolism, we kept our model at simple metabolic level. The proposed stoichiometric model accounts for the reactions of glycolysis, glyoxalase pathway, and thiol metabolism ( Table 3). It is well known that MGO is highly toxic entity and irreversibly modifies DNA, proteins/enzymes (Thornalley, 1993;Lo et al., 1994;Westwood et al., 1997;Greig et al., 2009) by reacting with Arg and Lys residues to form Advanced Glycation End products (AGEs) and promotes the formation of free radicals (Fraval and McBrien, 1980;Desai and Wu, 2008;Kalapos, 2008). Since, the formation of AGEs is a complex process in which MGO reacts with Lys and Arg residues of the proteins/enzymes, in our model this process was added in generalized form in which MGO directly reacts with Lys and Arg.
The exchange reactions were added to allow input-output exchange of the extracellular metabolites to enter the system from the medium or excretion of end products of any reaction out of the system. All reactions were set to the maximum for upper and lower bound limits for reversible/irreversible, uptake/excretion reactions except for glucose (Edwards and Palsson, 2000). For this small scale model, we took utmost care in including only the reactions those are required for maintaining the production of desired metabolites (Table S4). This model was submitted to BioModel Database (Chelliah et al., 2015) with the id: MODEL1909200001.

Formulation of Objective Function
We integrated the objective function (a mathematical equation providing drain of essential metabolites) consisting of the pseudo consumption of precursors from different pathways present in  The end product of glyoxalase pathway is D-lactate that ultimately leads to the formation of pyruvate. For simplicity, the stoichiometric coefficients of the precursors and metabolites used in formulation were kept as default. We considered, for our model, the end products of both the pathways, glycolysis (pyruvate) and glyoxalase pathway (D-lactate), contributing to the total biomass production. The system was investigated under two scenarios, first, functional GLO I, and second, nonfunctional (absent from the system) GLO I .

Scenario 1: Functional GLO I
It was noticed that glycolysis carried constant flux to produce pyruvate. However, glyoxalase pathway, comparatively, carried maximum flux for the production of D-lactate. It was also noticed that in presence of functional GLO I, advanced glycation end products (AGEs) formed at very minimal level. Although, both the pathways carried different flux levels, the maximization of OF was seemed to be majorly contributed through MGO synthesis and GLO I/II activities ( Figure 9A). Other metabolites synthesized at very minimal level probably due to smaller scale of our model and lack of other reactions for the compensation of their consumption (Chauhan and Singh, 2019).

Scenario 2: Non-functional GLO I
In, second scenario, we blocked the flux through glyoxalase pathway by setting GLO I flux to "zero" making it inactive. In this case, we observed that glucose entered the system at same rate as in scenario 1. Even the flux rate was found to be more or less similar for glycolysis but only upto the synthesis of GA3P and DHAP. In fact, the flux through TPI enzyme was more than that of scenario 1 that probably contributed to the synthesis of MGO. It is noticeable that MGO was produced at twice the rate than the previous case. This higher level of accumulation of MGO was found to contribute to equal level of the formation of AGEs and MGO •− free radicals ( Figure 9B). It was worth noticing that in case of inactive GLO I, the flux moved from glycolysis to MGO synthesis. Although in real network system, pyruvate production will be happening at normal rate except for that the inactivation of GLO I should lead to the increased level (∼1.5 fold) of accumulation of MGO .

Systems Pharmacology Modeling
Initially, the basal kinetic model (non-inhibited model) was constructed by incorporating the reactions from glycolysis, Glyoxalase system and T[SH] 2 synthesis via arginine, glycine, glutamate and methionine. Further, considering the important role of T[SH] 2, its demand reactions were also included in the form of TXN, TDPx, and TryP. As mentioned in earlier sections that due to the highly reactive nature of MGO, it irreversibly modifies Arg and Lys residues proteins/enzymes (Thornalley, 1993;Lo et al., 1994;Westwood et al., 1997) leading to formation of AGEs and free radicals (Fraval and McBrien, 1980;Desai and Wu, 2008;Kalapos, 2008). Hence, these reactions were included in generalized form in which MGO directly reacts with Lys and Arg, and AGEs were represented in the form of modified Arg and Lys (MGO-Arg and MGO-Lys, respectively). The model also consisted of various reactions for free radical metabolism. Further, to examine the effect of inhibited GLO I reaction on the system, a perturbation was imposed using previous studied competitive LmGLO I inhibitors, S-2,4-dinitrophenylglutathionylspermidine (DNPGS) and S-4-bromobenzylglutathionylspermidine (BBGS) (ki: 669 ± 57 and 0.536 ± 0.040 µM, respectively) (Ariza et al., 2006). Hence, two different perturbed kinetic models were constructed by introducing their inhibition reactions in the non-inhibited model. These inhibition kinetic models were submitted to BioModel database (MODEL1909200002 and MODEL1909200003, BBGS and DNPGS, respectively) (Chelliah et al., 2015). Analysis of the time course simulation of non-inhibited model has showed the maximum formation and degradation of MGO via glycolysis and GLOI enzyme, respectively. the maximum concentration of MGO was considered as initial concentration for both the inhibitors. GLOI reaction involves conversion of HTA to SDL. Although HTA is formed in reversible manner from MGO and T[SH] 2 non-enzymatically, we assumed that the introduced inhibition reaction, along with level of HTA, should also affect the level of MGO. However, when the comparison was made between non-inhibited and inhibited models, it was observed that even after the incorporation of inhibition reaction for GLOI, only HTA was accumulated in the system and not MGO (Figure 10, Figures S2, S3).
Comparison of metabolites' concentration and their production rates among the models revealed that the synthesis rate was higher in non-inhibited and DNPGS inhibited model as compared to BBGS inhibited model (Tables 4A,B). Both the inhibitors affected the level of few metabolites (free radicals) but BBGS was found to have greater effect. In BBGS inhibited model, the HTA accumulation was observed to increase by 1.6 fold in comparison to DNPGS inhibited model. Although, we did not observe any change in the level of Arg-MGO and Lys-MGO concentration in the system, however, the rate of their formation is increased by 1-3 folds in case of BBGS inhibited model. Through our kinetic modeling study, we are able to shed an important light on the mechanism of inhibition of GLO I reaction.
Further, to observe how GLO I inhibition has influence on the perturbed models; sensitivity analysis of the networks was carried out. For this purpose, GLO I reaction parameters, GLOI Vmax was used. Although perturbation did not greatly influence the accumulation of MGO in the system, the sensitivity analysis demonstrated that the inhibited GLOI Vmax affected other reactions fluxes. It seems like BBGS had more effect on other reaction fluxes than DNPGS when compared with the non-inhibited model. Between the two perturbed models, GLO I inhibition using BBGS negatively  affected GLO I reaction flux while DNPGS did not have any effect (Figure 11). Also, despite of no effect of inhibition on MGO level in the system, the rate of BBGS inhibited GLO I reaction influenced the flux of MGO formation from 2 to 10 folds.

DISCUSSION
To discover the robustness and reliability of a biological network it is important to understand its basic structural framework by quantitatively analyzing it through mathematical tools (Watts and Strogatz, 1998;Barabási and Oltvai, 2004;Dorogovtsev and Mendes, 2013;Weber et al., 2017a,b). This can help in uncovering the meaningful information by reducing the system to a smaller level to understand the relationship among its various components (Burago et al., 2001;Ollivier, 2009Ollivier, , 2011Ni et al., 2015;Sandhu et al., 2015). In current work, we motivate the use of systems biology driven approaches including Forman curvature and Forman-Ricci curvature measures along with FBA for the analysis of complex metabolic network of Leishmania parasite. Forman and Forman-Ricci curvature of master network has shown that nodes with high curvatures differ drastically in structure from the ones with low curvature values and appear to be in denser clusters acting as bridges between the subclusters in the network. This measure not only helped us to identify the relatively more important parts of the network but also provided the information on identifying the subnetworks. This gave us the idea to break the network in submodels to study the properties of various pathways at smaller scale and has helped us to detect important substructures inside the network belonging to particular classes of vertices or edges. Importantly, our analyses have pointed out several central nodes belonging to very common yet crucial metabolic pathways including glyoxalase pathway and T[SH] 2 metabolism pathways (Figures 5, 7, Figure S1). The significant importance of these high curvature nodes has also been demonstrated while performing the robustness analysis on the network (Figure 8). In fact, high correlation between Forman/Forman-Ricci curvature measures and common network properties has also demonstrated and pointed out the structural and functional importance of "M7" model comprising these metabolites ( Table 2). The assumptions made, here, are in congruence with the validated experimental literature and are an inherent feature of the model which signifies the importance of the metabolites undertaken to study the impact of the curvature measures. However, any modification of the connectivity matrices would certainly impact the curvature measures for a particular built up resampled model and would also impact the bridging of the metabolites.
It is worth noticing that several of these metabolites (MGO and T[SH] 2 ) have very important role in the regulation of redox homeostasis of the parasite. MGO is highly cytotoxic carbonyl and can irreversibly modifies DNA, proteins/enzymes (Thornalley, 1993;Lo et al., 1994;Westwood et al., 1997) causing the formation of advanced glycation end products (AGEs) and free radicals.
The neutralization of MGO is carried out by glyoxalase system comprising two key enzymes, GLO I and GLO II (Vickers et al., 2004;Greig et al., 2009;Wyllie and Fairlamb, 2011;Sousa Silva M et al., 2012). These two enzymes sequentially convert MGO into D-lactate in a two steps process: GLO I mediated isomerization of spontaneously formed MGO-T[SH] 2 conjugate, hemi-thioacetal (HTA), to S-D-lactoyltrypanothione, followed by GLOII mediated D-lactate formation (Thornalley, 1990). GLO I is reported as an essential and rate limiting enzyme in Leishmania. The essential role of GLOI has been demonstrated earlier in GLOI attenuated Leishmania donovani strain (Chauhan and Madhubala, 2009).
Our FBA result has also demonstrated similar information when we applied two different conditions to the model. Although, the OF we designed only serves as a tool to identify feasible flux across the model system, this has demonstrated the key role of GLO I. The non-functional GLO I has caused "zero" flux and led to several fold increased production of MGO (Figures 9A,B). The increased production of MGO in the system caused its accumulation due to under-determined flux in the downstream glyoxalase pathway and, hence, lack of its detoxification. The accumulation of MGO in the system consequently led to the formation of AGEs and MGO •− free radicals. Moreover, in scenario 2, glyoxalase pathway has also reflected the rate limiting step in very agreeable manner in the form of carrying the "zero" flux through the glyoxalase pathway reactions. These findings revealed the importance of GLO I and are in the agreement with the previous reports. Although, due to maintaining the simplicity of the model, we have not shown the effect of MGO •− free radical on other components in the system, however it has been demonstrated earlier that the excess amount of MGO causes the formation of increased level of AGEs and MGO •− free radicals (Nohara et al., 2002;Desai and Wu, 2008;Kalapos, 2008) and can inhibit the growth of Leishmania.
In our model, emergent metabolic properties of GLO(I) have been studied in a piecewise manner considering functional and non-functional GLO(I). We are well aware that there might be a change in metabolite vs. enzyme levels and any feedback introduced at all enzyme levels in the pathway of interest may give rise to different phenotype and correspondingly stable or sustained oscillations. Having said that, with our simplified model using a mathematical formalism wherein output depends on only two conditions in the pathway: functional and non-functional GLO I, these studies have been taken into consideration. Moreover, a measure of robustness in metabolic network accounts for the relative changes in the concentration of metabolites to all other metabolites. Our reduced model with local perturbations has a larger impact in a restricted subnetwork model system which is able to absorb the impact which a whole system model might not be able to bear.
Having stated above, for most part of the reconstructed network, taking into account the different point of observations, interaction dynamics does not drastically changed. The integrated analysis laid in this paper from computational perspectives helped pinpoint an interesting phenomenon that the underlined network model structure is largely preserved across the large common substructure models. Nonetheless, there might be instances at which major changes can occur in the network often in response to the metabolic or environmental stress. With the theoretical framework laid in this manuscript we are able to represent an in silico analog of the curvature measure and FBA. Therefore, the strategy outlined here may become an indispensable tool in future for analyzing smaller sub-network models created from a larger model. Using a model like this amounts to endowing the system with enhanced robustness and investigate how the dynamics might get affected. GLO I as a target deciphered in present work can further be validated through kinetic modeling and experimental validation to pinpoint its influence on direct or indirect metabolites.

CONCLUSION
In the quest to improve the current understanding of biological networks we have devised a newer strategy of formulation of curvature measure along with FBA in order to explain and predict biological targets. The methodology adopted in this paper has enabled us to determine the implicit relationship between metabolites. The model prototype developed in current paper largely depends on its structure and topological components. On the broader note, we have developed our mathematical model in an iterative fashion to optimize the flux of the metabolic system under limitation of certain enzymatic reactions thus leveraging the unifying framework for curvature measures along with FBA. GLOI deciphered as an important target both from curvature as well as FBA enabled our deeper understanding into the parasite redox metabolism.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this manuscript will be made available by the authors, without undue reservation, to any qualified researcher.

AUTHOR CONTRIBUTIONS
SS laid the design idea and implementation of the project. NC performed the computational modeling. NC and SS wrote the original manuscript and approved the final version of the manuscript.

ACKNOWLEDGMENTS
NC acknowledges her fellowship from DST-SERB NPDF (PDF/2016/001346) from Department of Science and Technology, Ministry of Science and Technology, Government of India. The authors would also like to thank the Director, National Centre for Cell Science (NCCS) for supporting the Bioinformatics and High Performance Computing Facility at NCCS, Pune, India. This manuscript has been released as a Pre-Print at bioRxiv (https://www.biorxiv.org/content/10.1101/ 512277v1; doi: 10.1101/512277.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fbioe. 2019.00336/full#supplementary-material Figure S1 | Depicting the internal connectivity in small scale models. The components are connected by backbone forming high curvature nodes (encircled in red color). Low curvature nodes are shown without encircled nodes. High curvature nodes act as bridges between major network communities and play important role in determining structural and functional organization of network.  Table S1 | List of various reactions and assigned rate laws used for kinetic modeling. Table S2 | Kinetic parameters and initial concentration of metabolites used in the kinetic model system.