The lower part of glycolysis is reconstituted in-vitro in a reaction assay medium described in a recent work (Moreno-Sánchez et al., 2008), containing different recombinant enzymes (PGAM, ENO and PPDK). The reaction was started by adding 3PG (4 mM). An additional reaction is added, the formation of lactate with lactate dehydrogenase (Figure 3A), in order to follow the flux of the overall pathway by following the rate of NADH oxidation, for more details, see Moreno-Sánchez et al. (2008) works. Concerning the peroxide detoxification pathway (Figure 3B), each enzyme was individually titrated, while keeping the other parameters in the in-vitro system constant. The pathway flux was determined in parallel by observing NADPH oxidation, see González-Chávez et al. (2015) for more information. Finally, the experimental procedures that were followed to obtain penicillin production data are described in the studies of Goldrick et al. (2015).

Two datasets are constructed here by applying data augmentation, using a gray-box model detailed in one of the following sections. For the first one, an exploration around the experimental data flux (43 ± 10 nmol·min⁻¹) from Moreno-Sánchez et al. (2008) at pH 6 is conducted. In fact, a sample of 2,000 normally distributed enzymatic balances was generated with the sample function on RStudio and resulted in a predicted flux between 0 and 53 nmol·min⁻¹ with the gray-box model. The term balance refers to a set of activities of the enzymes involved in the cascade of reactions. The second dataset is made up of experimental and predicted (gray-box model) data of PGAM, ENO and PPDK activities and pathway flux (J). The experimental data are obtained from plots of Moreno-Sánchez study (Moreno-Sánchez et al., 2008) (only the dots), while the predicted data are obtained with the gray-box model developed in a recent work (Lo-Thong et al., 2020), by varying each enzyme activity from 0 to 1000 mU with a step of 25 mU. These datasets are shown in Supplementary Tables 7, 8 respectively.

The second studied pathway consisted first of 58 experimental enzymatic balances and their corresponding flux. After applying data augmentation by using a gray-box model of this pathway, a bigger dataset of 1,671 data was obtained. As with the previous dataset, a combination of data normally distributed is generated with the sample function on RStudio, resulting in a predicted flux ranging from 0 to 11.46 nmol·min⁻¹. The new dataset is a mix of the previous experimental data and new predicted data of enzyme activities (TryR, TXN and TXNPx); final flux and is shown in Supplementary Table 9.

The two following pathways are modeled with an open-source software called COPASI (Version 4.24) (Hoops et al., 2006): the second part of glycolysis and the peroxide detoxification pathway. This software is used for metabolic network design, analysis and optimization. The first gray-box model, representing the lower part of glycolysis, is taken from Lo-Thong et al. (2020) work. It is based on the use of enzyme properties, including kinetic parameters and kinetic equations. To enhance the flux predictions, they suggested adding an adjustment term to the PPDK kinetic equation. The whole process concerning the composition of this term is explained in the previous work (see Methodology part of Lo-Thong et al., 2020 and Supplementary Table 1).

The second gray-box model represents the peroxide detoxification pathway and is built specifically for this study. It contains kinetic parameters and equations of three enzymes: TryR, TXN and TXNPx (Table 1). Also, we proposed to add two adjustment terms in TryR and TXNPx equations to improve flux predictions (Table 1). These are determined in the same way as the terms used for the glycolysis pathway. In fact, a first model was provided by González-Chávez et al. (2019) and could predict the final flux quite well when TryR and TXN activities were varied. However, it overestimated the flux when TryR activity was varied and underestimated it when TXNPx activity was varied. Therefore, we suggest adding a first adjustment term α(V_f − V_f0) in order to increase TryR rate and a second adjustment term β(V_f − V_f0) to decrease TXNPx rate. In these adjustment terms, α and β are defined numbers selected as the best for flux prediction from a tested range, V_f is TryR (or TXNPx) maximum rate in the forward direction in the model and V_f0 TryR (or TXNPx) maximum rate in the forward direction used in the in vitro reconstitution. Also, as V_f of TryR (or TXNPx) is equal to V_f0 when TXN's/TXNPx's (or TryR's/TXN's) activity is varied, we multiplied α (or β) by V_f − V_f0, so that the adjustment term would be zero when V_f = V_f0 and the flux predictions are not modified in these cases mentioned above.

Also, residual values are determined to evaluate how accurate the gray-box model is, and calculated as follows (1):

where e is the residual, y is the observed value and ŷ the corresponding predicted value.

For the datasets with <100 data, a process called data augmentation is performed. It consists of using models that accurately predict the experimental data to generate a new bigger dataset. Two different gray-box models are used in this study for the lower part of glycolysis pathway, retrieved from a recent study (Lo-Thong et al., 2020), and for the peroxide detoxification pathway (built for the present work). The gray-box models built on COPASI is set up to predict the variation of the final product concentration over the first hour for a given set of enzyme activities; then the COPASI outputs are processed to obtain the final flux of the studied metabolic pathway. Also, the overall process from the one-hour simulation for each enzymatic balances to the determination of the final flux is then automatized and applied to a range of enzymatic balances detailed in the previous subparts (Lower Part of Glycolysis Datasets and Peroxide Detoxification Datasets).

A brief analysis of the datasets is performed, including an examination of data distribution and the calculation of linear correlations between the input and output variables.

The determination of linear correlation between the inputs and output variables allows the assessment of the non-linearity for each studied metabolic pathway. As a rule of thumb, we consider that the non-linearity is high when one or more inputs has a linear correlation lower than 0.6. The lower the linear correlation, the greater the degree of non-linearity of the pathway.

To model the metabolic pathway, different machine learning models are developed on RStudio (Version 1.2.5001), with the help of Classification And Regression Training (caret, Version 6.0-86) (Kuhn, 2020).

The datasets are split into 80/20 for the training and test sets, and a k-fold cross-validation (with k = 10 for Dataset 1, 2 and k = 3 for Dataset 3) is performed on the models with the training set.

After this, the best models are selected based on:

The root-mean-square error (RMSE):

with Y_i and Ŷ_i being respectively the observed and predicted values, n being the total number of values and i = 1, 2…n;

the coefficient of determination (R²):

with Y_i and Ŷ_i respectively the observed and predicted values, n being the total number of values and i = 1, 2…n.

Also, a calculator was used for modeling the metabolic pathways, which has the following characteristics: cluster 2x Intel Xeon E5-2630v4 Broadwell-EP @ 2.20GHz 10 cores, 8x 16GB of RAM, 2400MHz, DDR4, ECC.

In addition, another type of metabolic pathway we can examine is the production pathways; their modeling would allow the development of an optimized overall process. In fact, another study revealed that ML methods can accelerate the optimization of chemical synthesis (Hein, 2021). As stated before, we do not need to enlarge this dataset, which is composed of records of the various parameters of an industrial-scale penicillin fermentation process. The use of this dataset made of only experimental data will ensure the reliability or not of the ML models for metabolic pathway prediction. It is important to consider that the inputs of our models are no longer the enzymatic activities, but different variables such as: batch time, oil flow, aeration rate, vessel volume and weight, carbon evolution rate and CO₂ percentage in off-gas. A slight variation of CO₂ in off-gas is recorded (Supplementary Table 9); this can be explained by the implementation of a system, by the operators, allowing corrective measures to be taken when the CO₂ level is too high, thus avoiding the detrimental effect of an accumulation of CO₂ on the growth of Penicillium chrysogenum and the production of penicillin. As the percentage of CO₂ in off-gas is maintained at a certain level, it is not surprising that the carbon evolution rate does not vary much either and presents a low standard deviation (Supplementary Table 5). Also, the output we are interested in is not the pathway flux, but the final concentration of penicillin (Figure 3C). As regards the correlation coefficient between the variables, we note that it is generally high between the parameters and the final penicillin concentration (Table 4); this correlation can be positive (e.g., time) or negative (e.g., oil flow). These correlation coefficients reveal the linear nature of the fermentation process studied in Dataset 3.

After showing that non-linear ML methods are more suitable for modeling metabolic pathways, we performed a comparison of the performance of all models. At first glance, the plots further confirm the preceding results and display higher RMSE values and lower R² values for the linear models compared to non-linear models (Figure 9). In addition, regardless of the number and/or type of data, we observe that Spike-and-slab, Ridge, Lasso and Bayesian GLM models give almost the same results (Figure 9 and Table 3). Also, it appears that some non-linear models work less well with large datasets; this is the case for ANN, bagEarth GCV, SVM Poly and SVM Radial (Figure 9). Moreover, it appears that random forest models (QRF and Rborist) are the best suited for metabolic pathway modeling, as they give the best results in term of RMSE and R² whatever dataset was used. Furthermore, we can evaluate the impact of the degree of non-linearity of the pathway on the predictions. Indeed, the pathway that has a high non-linear structure (Dataset 1) gives worse results for linear models than the pathway that presents a less non-linear structure (Dataset 3), which also gives good results with non-linear models (Figure 9A and Table 3). For example, Dataset 1 performs less well with the Ridge model, with RMSE = 29.311 nmol·min⁻¹ and R² = 0.716, than Dataset 3, which performs well with the same model, with RMSE = 3.579 nmol·min⁻¹ and R² = 0.87.

Besides, with a view to applying these methods at an industrial level, we perform a comparison of model error prediction and time of processing among the different datasets (Figure 10). The results confirm the previous findings, where random forest models have the best performance for metabolic pathway flux prediction. We noted that Rborist model presents a better RMSE - time of processing ratio than QRF model. However, even if QRF models have a processing time higher than 1h, we obtain an RMSE gain of about 96 %, when comparing it with PLS model, which could be of considerable significance for the industrial level. In view of the considerable gain of using this method compared to a linear one, non-linear methods could be more beneficial at the industrial level, where a gain of 1% is colossal. Spike-and-slab, Ridge, Lasso and Bayesian GLM models result in comparable performance in terms of RMSE and time of processing. At least, these results show a better RMSE – time of processing ratio for non-linear methods than for linear ones. We did not add the ANN models in the results, as they were not performed using parallelization process compared to the other methods.

Furthermore, we assess the impact of the amount of training data on ML model performance to have a desired level of performance (Supplementary Figure 9). We observe that the results are roughly the same for the datasets when they are predicted with linear models (Supplementary Figures 9B,D,E,H,I), thus the amount of data required to obtain a strong linear model can be higher than 80,000 data, as long as the studied pathway does not have a high degree of non-linearity. When it comes to non-linear models, we find that using a dataset smaller than 40,000 data is sufficient to obtain a good ML model (Supplementary Figures 9A,C,F,G,J–L). Using a dataset higher than 40,000 data leads to non-linear models that are efficient only in case of random forests (QRF and Rborist), Cubist and XGboost Linear methods, for which RMSE is low. We could also consider making an ablation of our datasets to examine the impact of amount of training set data on the ML model performance.

The first objective of this study is to determine what sort of data-driven model could better simulate the biological pathways studied. By using different datasets, we build several models with the enzyme balances or parameters collected from a bioreactor and reveal that Random Forests (QRF and Rborist), Cubist and XGBoost Linear are three good methods to predict the final flux or concentration of a final product. This works is part of a larger study about the applicability of either a knowledge-based or a data-driven approach. Indeed, in other fields such as fault detection and diagnosis, a comparison of these two methods demonstrates that they both have comparable performance and can be used (Alzghoul et al., 2014; Yang and Rizzoni, 2016). In biological system modeling, as is the case here, we demonstrated that in instances where little knowledge is available and difficult to obtain on a large scale basis (e.g., kinetic parameters k_cat and K_m of an enzyme, pathway fluxes), or when complex feedback regulation mechanisms take place, a data-driven method can be a good alternative for modeling a metabolic pathway, as many authors have shown before (Ramachandran et al., 2011; Hou et al., 2016). By comparison, the knowledge-based method can be laborious and long, due to data mining from the literature or wet laboratory experiments, whereas there is ease and speed of building models with the data-driven method (Kadarmideen, 2016).

Another criterion that we considered was the degree of non-linearity of the pathway. As mentioned above, it is generally admit that metabolic systems have an inherent non-linear behavior (Koza et al., 2001; Song and Ramkrishna, 2013; Yasemi and Jolicoeur, 2021). However, there is no formal demonstration of the non-linear structure of metabolic pathways. According to Song and Ramkrishna, this non-linear behavior would be due to: (i) the non-linearity of the chemical reactions forming the pathway and (ii) the regulatory processes that added non-linearity to the system (Song and Ramkrishna, 2013). Also, it is expected that pathway fluxes are non-linear, because they are controlled by enzymes and the activities of metabolic enzymes are saturable by their ligands. Besides, when the fluxes are measured in intact cells, they give a non-linear behavior and flux variation appears as hyperbolic or even sigmoidal. If the measured fluxes appear linear, it might be because the saturation point is not reached. Furthermore, according to the Metabolic Control Analysis, the fluxes are hyperbolic or non-linear because always exist one or two flux-controlling steps which utlimately determine the pathway flux (Fell, 1992). The determination of linear correlation coefficients of the different variables of the datasets gives us insights into the degree of non-linearity of the studied metabolic pathways and provides a method to evaluate the non-linearity of metabolic pathways. We found that all metabolic pathways studied here have a notable non-linear structure, with Dataset 1 having the highest degree of non-linearity, then Dataset 2 and lastly Dataset 3. These results generally comfort the main hypothesis that metabolic pathways are predominantly non-linear. The determination of the degree of non-linearity is therefore important for selecting and applying of a ML technique when modeling a metabolic pathway.

Moreover, the suitability of using either method relies on the quantity and quality of the knowledge or the data. Here, to illustrate this point, we simulate two datasets: the first one consisting of an exploration of the experimental data (2,000 data) and the second one composed of enzyme activities from 0 to 1,000 mU (68,950 data). The largest one gives better predictions for the three best models (Random forests, Cubist and XGBoost Linear) than the other dataset, and shows us the importance of having a large dataset before using machine learning methods. In fact, the size of the training set has been shown to be a major driving factor of prediction accuracy (Somarathna et al., 2017). However, we used two datasets made up of a mix of experimental and predicted data to build the models, and even if predicted from a good quality model, they remain mostly predicted data and are not comparable to a fully experimental dataset, which is also difficult to obtain. Thus, it would be worth considering methods using only experimental data, when sufficient data are available to build the models. Interestingly, a data-driven approach is often used to discover biological pathways or unravel pathways that are not well understood. Thus, combined with the knowledge-based approach, this can quickly make clear the complexity of biological systems modeling. Another possibility would be to test ML models on experimental data from E. coli or yeast, which can present a larger degree of non-linearity and are easily found in the literature. This issue will be addressed in our next study.

Surprisingly, model performance was weaker for the largest dataset from the bioreactor records than for the smaller datasets. The reason for this result may lie in the choice of input variables. Several studies have highlighted the need for variable selection in order to have better predictions (Camacho et al., 2018; Awan et al., 2019; Genuer et al.). Indeed, variable selection allows the use of the most informative variables to predict the output variable(s) and reduce the time of computing. Unlike the knowledge-based model, a diversity of variables for data-based models does not always mean better performance. This is one of the limitations of our study, since only one combination of input variables was tested during the work. It would be interesting for a future study to compare, for the same dataset, models using different sets of input variables, and to analyze their impact on model effectiveness.

Another major issue facing users of machine learning approaches is the interpretability of these models. Even if, at this time, we do not have a common general definition of this term, many researchers, such as Schmidt et al. (2019), define a model's interpretability based on two aspects: (a) intrinsic interpretability (or transparency): the ability to understand the inner mechanism of the model in the context of the study (e.g., identification of variables most involved in the predictions), and (b) post hoc interpretability: the ability to extract new information from the model or provide new insights into the relationships discovered during the process (e.g., the effect of a variable on another one) (Murdoch et al., 2019; Schmidt et al., 2019; Pintelas et al., 2020). Although some ML methods, such as decision trees or linear regression models, are easily interpretable; this is not the case for most of the models developed here (e.g., XGBoost Linear, bagging MARS, ANN). Nevertheless, using the variables that are directly related to the variable to be predicted, as we do here, allows us to gain some understanding of how the model works and the types of relationships that are revealed, enabling us to rely on the models. Furthermore, while we identified Random forests as one of the best methods for predicting final flux or product concentration, Pintelas et al. (2020) classifies it as a model that is hard to interpret. Therefore, it would be interesting to compute variable importance or to apply different techniques to explain the model in order to increase its interpretability (Zhou et al., 2019; Azodi et al., 2020). Besides, knowing that models based on decision trees are among the simplest to interpret, we support the idea of Schmidt et al. that RF models are more accessible than others from an interpretability point of view (Schmidt et al., 2019). An alternate solution would be to develop simpler models, but this would certainly reduce their overall performance.

Moreover, one of the key factors in the interpretability of the models is linked to the equations used. In fact, compared to knowledge-based models that use well-defined equations with a biological significance, ML models are governed by other equations, which sometimes are “outside our understanding” as Schmidt et al. (2019) observed in their study of the applications of ML in solid-state materials science. This raises a real problem of confidence in the prediction results obtained with such methods. As these authors point out, the fact that these models were not based on physical principles in their studies, or on biological principles in ours, could result in wrong predictions in completely unexpected cases, while providing great results overall. And in the present case where the models are used in the context of biomarker identification or optimization of an industrial bioreactor, we cannot risk obtaining such results from our models in these specific situations. Far from hindering us in the use of ML models, awareness of these problems allows us to formulate several recommendations for future research. These include the combination of interpretable models, e.g., knowledge-based kinetic models with ML models, e.g., random forests models; the prediction of a new set of experimental data with unexpected values. In this latter instance, this would require experimentally testing a range of “extreme” data that would be found in the parasites studied, or recording the bioreactor data even during failures of the penicillin production.

After analyzing the interpretability of the different modeling methods, it is worthwhile to note some advantages and disadvantages of their use in flux and concentration prediction. One of the best methods in our case is the random forest (QRF and Rborist). Many studies report the use of random forest in the biological field for the prediction of: protein interaction (Qi et al., 2006), drug response based on protein markers (Ma et al., 2006) and in vitro drug sensitivity (Riddick et al., 2011). Also, Riddick et al. used SVM and random forest to predict the flux of N₂O emissions, and found that random forest achieves the best performances among the built models (Villa-Vialaneix et al., 2010). They highlighted that these models offered the advantage of having a low computational cost, compared to the SVM method. However, in our case, we notice that random forest is the least accurate predictability model compared to SVM methods, with the highest computation time for almost all datasets. Moreover, among the random forest packages developed on R, Rborist is quite a recent implementation, designed for multicore hardware, which minimizes data movement within memory to increase the performance and decrease the processing time (Wright and Ziegler, 2017). Surprisingly, here, Rborist package is the one that has the longest time of computation and is more efficient on big datasets compared to other methods. It would be of interest to create variant models combining the random forest method and other methods, as in previous studies (Chen et al., 2018; Zampieri et al., 2019). An existing variant of random forests is the quantile regression forest (QRF) method, which has the capability of establishing prediction intervals that cover uncertainties, useful in the prediction of possible new data (Meinshausen, 2006). Francke et al. demonstrated in their work that this method had the advantage of calculating uncertainties associated with the predicted sediment yields, through the calculation of confidence intervals (Francke et al., 2008). But they also stated that the model predictions will always be within the range of observations, which prevents implausible values but inhibits prediction outside the range of values learned from the training set. We saw here that, overall, QRF models have a good generalization capability; additional prediction of new experimental data, with data separated by a larger stepsize (>25), would be beneficial to confirm or invalidate this capability. This could be useful for the study of metabolic pathways in extremotolerant organisms.

This leads us to note one of the advantages not only of the QRF method but also of other ensemble learning methods, such as XGBoost Linear: prediction from high-dimensional data. Indeed, these models are among the best we have, with any starting dataset we have, from the simplest to the most complex with several types of variables. Remarkably, compared to other models, XGBoost Linear is better ranked for small datasets. This is confirmed by the work of Yang et al. (2010) which propose that ensemble methods have the advantage of reducing the potential for overfitting in small sample size problem. Another strength of XGBoost Linear compared to its peers is the combination of high accuracy and a short time of processing. However, despite the great accuracy of these models, they are often more complex and less interpretable, and present a higher computational intensity.

Moreover, Cubist, a model based on modified regression tree theory, has the advantage of analyzing big data with high speed (Xu et al., 2018). This was confirmed by our results, which show that Cubist is one of our best models (e.g., for Dataset 1, Cubist: 2.49 min and QRF: 1.76 hr). However, we noted that the performance was better for the small datasets than for the bigger one. Another advantage that Das et al. noticed is the fact that the Cubist model is easy to interpret and is a suitable method for beginners (Zhou et al., 2019; Das et al., 2020).

The PLS method turned out not to be appropriate here to model these pathways and predict the final flux starting from enzyme activities, or the final product concentration starting from parameters of a bioreactor. This may be due to the inherent limitation of the PLS method to capture the non-linearities of the metabolic pathways. However, it performs better when we have a smaller dataset, as it has also been noted in a previous study on gluconeogenic flux prediction (Antoniewicz et al., 2006). But these results contradict those obtained with the PLS model for the prediction of limonene and isopentenol synthesis. In fact, in this work, results showed that the model performed well when the dataset was larger (lower RMSE, better predictions) (Costello and Martin, 2018). Also, one big advantage of the PLS technique remains that it has the shortest calculation time for modeling.

It is relevant to observe that the model implementation will differ depending on varying levels of data. In fact, a ML model will be more difficult to implement, if the available data is limited. In this case, a significant additional time is required. Among the various studied models, the difficulty to implement the model could also be based on the higher number of parameters to adjust during the training time.

Our findings generally support the idea that non-linear models are more suitable than linear ones for modeling metabolic pathways. Furthermore, it would be interesting to apply these ML models on genome-scale metabolic networks for which the literature abounds in data. Recently, hybrid models coupling a genome-scale model and ML model have been found to be effective for different purposes such as the prediction of individual amino acid concentration in culture medium (Schinn et al., 2021) and identification of prognostic metabolic biomarkers in cancer studies (Lewis and Kemp, 2021). One of the benefits that ML models could bring is the integration of multi-omics data as genomic, transcriptomic, metabolomic and proteomic data. This topic will be addressed in an upcoming study.

As far as we know, genome-scale models have a predominant place in the field of metabolic networks for the identification of key-molecules in the metabolism. This study allows us to consider the machine learning methods as performant models to predict metabolic pathways. Indeed, their ability to take over large datasets makes them applicable techniques to efficiently predict larger metabolic pathways (e.g., E. coli). While flux balance analysis (FBA) based methods, as used in the genome-scale models, need information about the pathway in a given condition as they are hypothesis-driven, machine learning models could predict the metabolic pathways without needing to clearly understand the underlying biological mechanisms of the pathways. Also, constraint-based model (e.g., FBA) are not able to predict metabolite concentrations, while the machine learning methods can consider these predictions. We can thus envisage a hybrid method using both machine learning and FBA methods for metabolic pathway modeling (Zampieri et al., 2019).

Given the many different methods that exist and continue to emerge, one can struggle with the choice of a model to build from a dataset. Faced with this decision, we can choose to build simple models or to use models being used in the same field of study and giving good results (Camacho et al., 2018; Cifuentes et al., 2020). In view of this, it would be useful to review and define some basic rules for building a decision-making support for future studies on modeling metabolic pathways. The first feature to consider is the quality of the biological dataset (Figure 11). Do we have an initial dataset of good quality? Data quality can highly impact the model predictions. If the model is not of good quality, it would be better to build a new dataset and generate good quality experimental data. When the dataset is of good quality but small in size, it is useful to do data augmentation, as we did in this study; if this is not possible, we can use an ensemble model to build the metabolic pathway, since such models can deal with small datasets. Another useful criterion we can investigate is the number of variables. If the dataset presents a high number of variables, we can consider doing variable selection before building the model, or we have the option of building the model by using ensemble modeling that gives good accuracy with several input variables. Also, one key factor is the non-linearity of the studied metabolic pathway; do we have a non-linear or a linear process? If our pathway is linear, we can design a battery of linear models which will give a high performance. But if our study involves a pathway that is non-linear, then it is preferable to use a non-linear model. After building our model, an evaluation of its accuracy is necessary to validate it. In case the performance of the model is not suitable, we can plan to refine it, for example by tuning the hyperparameters (Chicco, 2017), or simply to replace it and build a new one.

Non-linear machine learning methods enable us to model metabolic pathways by identifying key-molecules, which are important for the drug-design process, improving disease diagnosis (cancer, viral/parasitic/bacterial infections, neurodegenerative diseases) by highlighting the differences between healthy and pathological situations, or even optimizing industrial production processes.