Hybrid Dynamic Models of Bioprocesses Based on Elementary Flux Modes and Multilayer Perceptrons

: The derivation of minimal bioreaction models is of primary importance to develop monitoring and control strategies of cell/microorganism culture production. These minimal bioreaction models can be obtained based on the selection of a basis of elementary ﬂux modes (EFMs) using an algorithm starting from a relatively large set of EFMs and progressively reducing their numbers based on geometric and least-squares residual criteria. The reaction rates associated with the selected EFMs usually have complex features resulting from the combination of different activation, inhibition and saturation effects from several culture species. Multilayer perceptrons (MLPs) are used in order to undertake the representation of these rates, resulting in a hybrid dynamic model combining the mass-balance equations provided by the EFMs to the rate equations described by the MLPs. To further reduce the number of kinetic parameters of the model, pruning algorithms for the MLPs are also considered. The whole procedure ends up with reduced-order macroscopic models that show promising prediction results, as illustrated with data of perfusion cultures of hybridoma cell line HB-58.


Introduction
Mammalian cell cultures are now widely established to produce therapeutic products of interest such as monoclonal antibodies, viral vaccines and proteins used in the treatment of genetic diseases.Following the current trends of Process Analytical Technologies (PAT) and Industry 4.0, the development of dynamic mathematical models and predictors or estimators is receiving ever-increasing attention in view of establishing digital twins of the bioprocesses.
The models can usually be classified into structured and unstructured models, where the first group considers the cellular metabolism, while the second neglects the intracellular activity and focuses on the evolution of extracellular species.Over the years, connections have been established between the two approaches, and intra-and extracellular information can be combined to develop reduced-order macroscopic models [1][2][3][4].A central concept in the development of minimal bioreaction models is the notion of elementary flux modes (EFMs) [5], which can be seen as the simplest metabolic pathways linking extracellular substrates to products.The EFMs form a convex basis of the flux space, whose dimension unfortunately grows very quickly with the size of the metabolic network.Tackling this combinatorial explosion has been a major research concern in the last decades.
To alleviate this issue, a first approach is to consider metabolic networks of modest sizes, which can be obtained on the basis of detailed networks through metabolic flux analysis [6], for instance, by discarding insignificant fluxes [7].In this way, elementary vectors can all be computed and enumerated using specific software tools such as METATOOL [8] and EFMTOOL [9].
Alternatively, when the number of EFMs is prohibitively large [10], methods to select EFM subsets are required.To this extent, Figueiredo et al. [11] identified the shortest elementary modes, Kaleta and his group [12] proceeded with subsystem analysis, Jungers [13] decomposed the flux distribution into a minimal number of elementary vectors, Machado and Tabe [14,15] developed random sampling, Soons [16,17] used ranking or controlled random search to select a subset of modes based on an optimization criterion and Oddsdottir [18,19] employed column generation techniques to identify a small number of modes.Most of the previous methods account for constraints related to cell-specific uptake or secretion rates to further constrict the flux space, as studied in [20,21].
Another strategy consists in reducing the number of elementary vectors from a larger initial pool so as to keep only a small set of representative modes.In connection with this, Hebing et al. [22] used an EFM reduction procedure based on a geometrical collinearity criterion.More recently, several procedures of EFM reduction have been developed by our research group, i.e., Abbate [23] selected the best EFM candidates based on the formulation of a linear optimization problem and Maton [24] developed a reduction methodology based on a combination of several criteria based on collinearity and a series of constrained least-square problems.
What makes the use of elementary flux modes attractive in bioprocess modeling is the possibility to derive minimal bioreaction schemes.Henceforth, it remains to model the kinetics using specific functions (such as Monod or Haldane laws or more complex nonlinear functions of the extracellular species) and parameter identification to establish a functional model from sets of experimental data.In the last decades, many studies addressed the issue of kinetic identification.However, the determination of kinetic structures is often based on arbitrary choices.Indeed, specific kinetic phenomena such as activation, saturation and inhibition can sometimes be described by different mathematical expressions with similar evolutions with respect to the culture species.Attempts have been made and achieved to propose general kinetic formalisms, as in [25,26], where power-law equations were used to represent overall reaction rates of complex biological systems [27][28][29][30][31][32].Nevertheless, although these methods are effective and have the merit of being systematic, they do not capture a double component effect.More recently, generalized Monod equations were proposed in [20,[33][34][35]; a systematic procedure to select the most likely kinetic structures using decision graphs, nested models and likelihood ratio tests was developed in [36]; and another systematic procedure to model the kinetics using generalized-kinetic functions and a three-step parameter identification procedure was detailed in [37,38].
All the previous approaches are based on kinetic functions incorporating factors representing various phenomena of activation, inhibition and saturation, and providing some biological interpretation.When the reaction rates are difficult to describe as they include influences from several culture species, an alternative approach is to resort to purely data-driven techniques such as neural networks [39][40][41].Although these models present a high flexibility to capture nonlinear kinetic phenomena, they lack biological interpretation.This latter situation corresponds to the starting point of this study.Indeed, the application of our EFM reduction algorithm [24,42] yields, on the one hand, a candidate stoichiometric matrix and, on the other hand, the time-evolution of the metabolic fluxes corresponding to the selected EFMs.These fluxes may have a complex time evolution resulting from the combination of several phenomena triggered by different culture species.Therefore, the objective of the present study is to explore the potential of recurrent neural networks, and particularly multilayer perceptrons (MLPs), to describe these fluxes.The main reason for this choice is their simplicity and the existence of efficient toolboxes such as the Deep Learning toolbox in Matlab.A disadvantage of MLPs, however, is that they have a fully connected structure with a quick increase in the number of parameters with the number of neurons in the hidden layer.In this study, it is desired to keep the number of parameters at the minimum and the use of a pruning algorithm is also explored.To support our findings, experimental data of perfusion cultures of hybridoma cell line HB-58 are used to derive a dynamic hybrid bioreaction model.This latter model is corroborated with the mechanistic model of [43], in particular, in order to interpret how overflow metabolism is represented by the neural structure.
The paper is organized as follows.The next section presents the concept of hybrid modeling using elementary flux vectors for the derivation of minimal bioreaction models and neural networks for the description of the kinetics.In Section 3, information is given about cultures of hybridoma cell line HB-58, which are used as a representative case study.Numerical results are discussed and interpreted in Section 4. Finally, conclusions are drawn in Section 5.

Hybrid Modeling
The hybrid modeling methodology is sketched in Figure 1.Starting from the definition of a metabolic network with stoichiometric matrix N, the measurement configuration N e and experimental uptake and excretion rates ν m , this method first exploits mechanistic information based on the concept of elementary flux modes and an EFM reduction procedure to infer a minimal bioreaction model and its corresponding stoichiometric matrix K. Second, the procedure exploits the time evolution of the reaction rates Φ num , which are provided by the EFM reduction algorithm, and applies a data-driven modeling of these rates using MLPs with the concentration of extracellular measurements ξ as inputs.The model components K and Φ can then be used in a predictive model.The full methodology is described in detail in the sequel.N is the stoichiometric matrix, N e is the stoichiometric matrix of extracellular measurements, ν m is the vector of experimental uptake and excretion rates, K describes the stoichiometry of the bioreactions, φ num is the vector of reaction rates obtained numerically from the reduction procedure, φ is the vector of reaction rates modeled using NN and ξ is the vector of concentrations of culture species.

EFM Selection 2.1.1. Metabolic Network Analysis
Cellular metabolism is defined as a set of chemical reactions, possibly catalyzed by enzymes, taking place within the cell and forming metabolic pathways.These intracellular reactions may be translated into a matrix representation defining a metabolic network as a m × n stoichiometric matrix, denoted N. In this formalism, m is the number of internal metabolites and n represents the number of reactions.Considering the pseudo-steady state assumption, the following homogeneous system of linear equations is obtained: where v ∈ R n gathers the fluxes of the network.Moreover, to express that the reactions have a net direction, the following constraint may be stated: Henceforth, the set of possible flux distributions v defines a pointed polyhedral cone S in the positive orthant.The edges of S represent elementary flux vectors-the simplest metabolic pathways connecting extracellular substrates to final products without accumulation of metabolites.
As depicted in Figure 2, it is possible to further reduce the solution space by adding a set of linear constraints making use of experimental uptake and excretion rates ν m : Subject to the constraint in Equation ( 2), the set of solutions is included in the polytope F (F ⊂ S) and only specific combinations of modes provide a solution.In Equation (3), N e is the stoichiometric matrix of extracellular measurements, a m e × n matrix where m e stands for the number of measured extracellular species.The vector ν m is inferred from the measurement of the time evolution of the extracellular concentrations, for instance, using regression based on the smoothed derivative of the concentration signals.

EFM Reduction Procedure
The EFM reduction method has been originally developed in [24,42], and only a brief overview is provided in this section.The 4-step method, illustrated in Figure 3, is divided into (i) the initial generation of the modes, (ii) the biological interpretation of the reaction scheme, (iii) a preliminary reduction up to Ω modes and (iv) the selection of a minimal set of Λ EFMs.

•
Generation of the initial EFM set: The first step concerns the generation of elementary flux vectors.If the size of the metabolic network is not too big, the whole set of modes can be computed and enumerated using software tools such as EFMTOOL.
Otherwise, if the number of EFMs becomes prohibitive, alternative methods to identify only subsets are required.A fast generation algorithm [13] can be used that requires also the knowledge of experimental measurements of uptake and excretion rates.Nevertheless, regardless the EFMs generation method, a matrix of elementary flux modes E can be obtained.

•
Biological interpretation: The second step consists in ensuring a biological interpretation of the matrix K, which is defined by and is a m e × n EFM matrix whose columns provide the stoichiometry of each bioreaction.All the modes conducting to a macroreaction with no reactant or no product should be removed, yielding a reduced matrix of modes, denoted E * .Note that reactions corresponding to maintenance are kept.

•
Main EFM reduction: This step allows reducing the number of modes up to a target Ω.This value is generally set close to, and sometimes slightly greater than, m e to avoid computational issues during the final step of the procedure.First, collinearity between vectors is evaluated and the collinear modes are discarded.Second, an optimization-based reduction is achieved where a randomly selected vector is removed if the following inequality is satisfied: where Ξ * is the performance index for the candidate elimination, Ξ is the prior value of the indicator and tol is a tolerance value defined by the user.The selected performance index is a weighted constrained least-squares problem: where ν m,k is the vector of uptake and excretion rates at every time step k, φ k corresponds to a time-varying decomposition of the flux ν m,k into a reduced set of vectors stored in K e , and W is a weighting diagonal matrix whose diagonal elements are max ).This criterion represents how well the positivity constraints can be satisfied while reproducing the evolution of the extracellular fluxes.Note that K e = N e E e , where E e denotes a reduced matrix of EFMs.

•
Selection of a minimal bioreaction model: For this final step, an even smaller number of EFMs is selected among the previous set of Ω modes (Λ < Ω).This target Λ is chosen below the number of extracellular measured species m e in order to derive macroscopic models with less reactions than components.This step is no longer based on random successive eliminations of modes but is a selection step of the best combination of Λ EFMs among the previous Ω modes.Hence, the performance index Ξ of all possible EFM combinations is computed and the final set of modes is the one with the smallest value of the indicator, which represents the distance to the experimental data.Note that the number of possible combinations is given by As a consequence, Ω has to be chosen small enough to avoid an unmanageable number of EFM combinations and, at the same time, large enough to have more flexibility in the selection, i.e., to have a sufficient pool of EFMs for the selection of the best EFM combination.

Dynamic Mass Balance Model
From the reduction procedure presented in the previous section, a stoichiometric matrix K, as well as the time evolution of the vector of reaction rates Φ(t) have been obtained.The product of these two quantities defines the biological uptake and excretion rates of the culture species ξ in the following mass balance equation system: In this equation, D is a diagonal dilution matrix and ξ in denotes the inflow concentrations.
To develop a predictive model, it is now necessary to represent the flux vector Φ by a parametric model describing the influence of the extracellular culture species.In this study, MLPs will be used to describe the kinetic laws.This weighted sum enters an activation function f to predict the output signal y = f (z).There exist many activation functions, with sigmoid functions being a popular choice.To represent complex nonlinear functions, the use of multilayer perceptrons (MLPs) is recommended, which are organized as follows: • one input layer, which distributes the input values to the first hidden layer; • one or several hidden layers of perceptrons; • one output layer, which recovers the output of each perceptron of the last hidden layer.
Figure 5 shows a typical multilayer perceptron with one hidden-layer in the context of the present study, i.e., the input signals are the species concentrations ξ i and the outputs are the elements of the flux vector Φ.It is worth noting that the number of hidden layers constitutes the true computational engine of the MLP and the more hidden layers, the more powerful the artificial neural network but at the expense of a large number of parameters and the risk of overparametrization.Nevertheless, even one hidden-layer MLP can approximate the mapping of any continuous function [44].The multilayer perceptron is known as a feedforward network as data flows in the forward direction from the input to the output layer.The neurons can be trained with a backpropagation learning algorithm computing the gradient of a loss function.Basically, the perceptron learning consists in adjusting the weights and the biases of the network in order to minimize the deviations between the outputs predicted by the network and the expected ones.For the training of MLPs, backpropagation is an efficient algorithm to estimate the gradients, which can be exploited in various methods, such as gradient descent or stochastic gradient descent, updating the parameters of the structure to minimize the loss function.More information on the learning algorithms can be found, for instance, in [45].
Before training the neural network, it is common practice to perform pre-and postprocessing steps on the network inputs and outputs in order to make the training process more efficient.Depending on the magnitude of the input values, sigmoid activation functions may become saturated leading to small gradients and slow convergence of the training procedure.To tackle this issue, a normalization step is applied to input signals in the data set during a pre-processing step and the output signals can be reversely transformed back into the original range in a post-processing step.
MLPs are fully interconnected and the number of parameters rapidly increases with the number of hidden layers and neurons in these layers.Pruning algorithms discard redundant and unnecessary connections [46][47][48][49][50].By removing synaptic connections with little influence, a sparsely connected network with multiple zero weights can be obtained.Methods to perform pruning of MLPs are proposed in [51].Figure 6 illustrates the architecture of a network before and after different techniques of pruning.In this study, a simplified version of the magnitude pruning algorithm proposed in [48] is developed.The idea is to assign a score to each parameter of the network (weights and biases) corresponding to its absolute value.Due to the pre-and post-processing, the score of each parameter is assumed to represent its relative importance to the accuracy of the trained network.Hence, a synaptic connection or a bias may be removed if its score is smaller than a threshold value defined by the user.

Case Study: Perfusion Cultures of Hybridoma Cell Line HB-58
Data used in this study come from the hybridoma cell line HB-58 enabling the production of IgG1 monoclonal antibodies, specific for mouse kappa light chain.Serum-free medium-prepared from DMEM/F12 (1:1) and completed with glutamine, 500 µmol of ethanolamine, 10 mg of bovin insulin, cholesterol, Pluronic F-68, HEPES and other additives-was used.The cells were kept as suspension cultures in the medium in shake flasks at 37 °C in a 5% CO 2 incubator.Experiments were conducted at the State Key Laboratory of Bioreactor Engineering, East China University of Science and Technology (ECUST), in Shanghai [52].
Perfusion cultures were performed in a 2 L stirred bioreactor and were settled in a working volume of 1.8 L during the whole duration of the culture.Cultures were carried out in a controlled environment (36.8 °C, 40% DO, pH = 7.0 ± 0.2).The perfusion phase started after 44-56 h of batch culture with a constant dilution rate of 0.0197 h −1 .Cells were retained by a spin-filter (20 µm) and the stirring speed was fixed at 200 rpm.Data acquisition and process control were performed using the supervisory software MFCS/Win 3.0.Table 1 summarizes the operating conditions of the cultures.

Metabolic Network
Depending on the level of description required, networks of different sizes might be considered.Henceforth, a metabolic network of hybridoma cells with 70 biochemical reactions and 44 internal metabolites is considered.The details of this metabolic network can be found in [53]; it includes the major reactions of central metabolism such as the pathways of glycolysis, the tricarboxylic acid cycle, amino acids metabolism, and the synthesis of biomass and antibodies.The pentose phosphate pathway is not taken into account in the network definition because, in most tissues, 80 to 90% of glucose oxidation is achieved by glycolysis.The stoichiometric coefficients of the biomass and antibody synthesis are taken from literature [52].Note that considering a metabolic reaction for the formation of biomass allows its prediction when deriving dynamical models.

Measurement Configuration
Only 6 extracellular measurements are accounted for-namely, glucose, lactate, glutamine, ammonia, alanine and biomass-gathered in vector ξ.This leads to the stoichiometric matrix of extracellular measured species N e .Furthermore, the experimental uptake and excretion rates ν m (t) may be computed using smoothing splines and differentiation methods.

EFM Selection
The first step of the reduction procedure consists in computing an initial set of elementary flux vectors.As stated in Section 2.1.2,different approaches exist depending on the size of the metabolic network.In this case, using EFMTOOL as the EFM generator, an initial set of 22,563 modes is obtained.Alternative methods such as the one proposed in [13] identify only representative subsets of EFMs and allow getting a few hundred modes as a starting point.However, in order to prove the effectiveness of the reduction procedure, the whole set of elementary vectors is considered in this study since their number remains manageable.Then, for the purpose of ensuring a biological interpretation of the reaction scheme, any vector leading to a macroreaction with no biological meaning is discarded and the reduced matrix of modes E * is made up of 21,874 vectors.The power of the methodology lies in the main reduction algorithm.With a collinearity indicator of 99% between the vectors, a drastic reduction to 214 modes is achieved.Next, an optimization-based reduction is performed with a target number Ω equal to 10.The last step is the selection of a few modes from the previous set to derive models with fewer reactions than components.In this study, only five elementary flux modes are finally selected (Λ = 5).Note that the target Λ could have been set to an even smaller value.However, because overflow metabolism is noticed in hybridoma cells, a further reduction may cause a loss of essential biological information.
A first direct validation of the reduced model can be achieved by comparing the experimental uptake and excretion rates ν m (t) to the product K e Φ(t) computed in the optimization problem.This validation can be pursued by analyzing the prediction of the measured concentrations by integration of the computed fluxes and identification of the most likely initial conditions of the measured species.These results are shown in Figures 7 and 8, respectively.Almost no loss of information is noticed and the prediction of the concentrations is very satisfactory, highlighting the merits of the modular EFM selection procedure.Besides the previous validations, it might be interesting to examine the final set of macroreactions obtained from the reduction algorithm: where α i , β i , γ i , δ i , i and σ i are stoichiometric coefficients, listed in Table 2. K is a m e × Λ matrix with m e = 6 and Λ = 5.Equivalently, a macroreaction scheme can be drawn: This reaction scheme makes sense and is validated by the study conducted in [43].Indeed, glucose and glutamine are consumed with biomass growth and metabolites production defining the so-called respiratory metabolism.Moreover, the phenomenon of overflow metabolism is pointed out with lactate production from glucose excess and production of ammonia, lactate and alanine from glutamine excess.This point will be further discussed in Section 4.3.

Kinetic Modeling
For the derivation of a dynamic model, it remains to describe the numerical signals Φ coming from the reduction algorithm by neural networks.As illustrated in Figures 1 and 5, a multilayer perceptron will use as input signals the concentration of extracellular species ξ and as target values the reaction rates Φ computed numerically in the reduction procedure.Therefore, there are six input signals, namely, the concentrations ξ i (i ∈ [1,6]) of the measured species and five output signals Φ i (i ∈ [1,5]), which are the specific rates of each reaction.Only one hidden layer is considered to limit the number of parameters (weights and biases), and the nonlinear activation functions in the hidden layer are hyperbolic tangent sigmoid functions while the activation functions in the output layer are identity functions.In summary, the following equations can be written: with This function has the same shape as the classical sigmoid function but horizontal asymptotes are in −1 and 1 so that the output values are squeezed into [−1, 1].In the previous equations, Z is the number of neurons in the hidden layer.Hence, the number of parameters is 5 + 12Z.Depending on the complexity of the problem, the number of neurons in the hidden layer can be increased at the expense of a significant increase in the number of parameters to identify.Table 3 summarizes the number of parameters of the network according to the number of neurons in the hidden layer.The rapid increase in the number of parameters with the number of neurons justifies the use of pruning algorithms.As a direct validation, Figure 9 shows the time evolution of the specific reaction rates Φ produced by the reduction procedure (and used as target values) and the prediction with MLPs with different numbers of neurons in the hidden layer.As expected, the more neurons, the better the fitting to numerical results.Figure 10 depicts the corresponding concentration profiles.Although significant deviations appear in the reaction rates, they have a limited impact on the prediction of the time evolution of the concentrations.Hence, an MLP with only three neurons in the hidden layer is an acceptable compromise between the fitting to experimental data and the number of network parameters.Note that experimental data has been partitioned to avoid overfitting, i.e., 70% for training, 15% for validation and 15% for testing.
In order to improve the MLP prediction, instead of backpropagation, a global nonlinear identification of the weights and biases of the network can be performed using the Nelder-Mead simplex optimization algorithm in order to minimize a weighted nonlinear leastsquares criterion: where Θ is the vector of parameters to be re-identified containing the weights and biases of the neural network, Φ ij gives the network prediction at the ith time instant in the jth experiment, Φ num ij is the vector of the corresponding numerical signals coming from the EFM reduction procedure and Σ is a normalization matrix in which diagonal elements are chosen as the squares of the maximum rate of each reaction.To improve convergence, the optimization algorithm can be run several times using the parameters found in the previous round as initialization.The confidence intervals and variation coefficients of the estimated parameters can be obtained from the inverse Fisher information matrix [54].The variation coefficients of the estimated parameters are given in Table 4.
Next, the magnitude pruning algorithm is exploited to reduce the number of parameters.The results of the pruning algorithm are shown in Table 5. W 1 is a matrix collecting the synaptic connections between the input layer and the hidden layer, b 1 represents the biases of the neurons in the hidden layer, W 2 contains the weights between the hidden layer and the output layer and b 2 are the biases of the neurons in the output layer.In this way, the number of parameters is reduced by 17% (41 → 34 parameters).
The time evolution of the reaction rates and the prediction of the measured concentrations using a pruned multilayer perceptron are illustrated in Figures 11 and 12, respectively.Figure 11 shows the added value of a global nonlinear re-identification of the parameters.Small deviations may be observed when using the pruned MLP for kinetic identification but the results remain very satisfactory, highlighting the merit of the simplified pruning method.Figure 12 validates the whole procedure providing good results for the prediction of the measured concentrations.
Furthermore, a pruning of the neural network can also be achieved on the basis of the coefficients of variation of the estimated parameters.Indeed, weights and biases identified with a coefficient of variation greater than 30% can be discarded without affecting the fitting of the experimental data.However, the elimination of the weights in W 2 must be exercised with more care and only parameters identified with a coefficient of variation greater than ∼ = 50% are neglected.Figure 13 shows the prediction of the time evolution of the extracellular measured species using this pruning strategy.It yields fewer parameters in the neural structure (41 → 30 parameters) at the expense of larger deviations from experimental data.However, the results remain acceptable for deriving a macroscopic model for control and optimization purposes.Furthermore, the pruning strategy using the value of variation coefficients covers up most of the parameters defined as nonessential when using the magnitude pruning algorithm, validating the two pruning approaches.

Interpretation of Overflow Metabolism
Overflow metabolism is a phenomenon in which cells achieve incomplete oxidation of an abundantly supplied energy source, commonly glucose and glutamine, despite aerobic conditions.It results in the excretion of organic end-products that are often inhibitory.Overflow of glucose leads to the formation of lactate while overflow of glutamine results in the formation of ammonia.As mentioned in [55], glutamine is mainly excreted as alanine, proline and aspartate.This phenomenon is essentially described by the reaction scheme deduced from the modular EFM selection procedure in Section 4.1.Indeed, in mammalian cell cultures, two main metabolic states can be pointed out: (i) the state of respiratory metabolism at low substrate uptake rates and (ii) the state of overflow metabolism at high substrate uptake rates.During respiratory metabolism, the consumption of glucose and glutamine leads to biomass growth and metabolite production such as lactate and ammonia.Respiratory growth is captured by reaction (14) for the substrate consumption and biomass growth, together with Equations ( 11), ( 13) and ( 15) for the metabolite production.Overflow of glucose leads to the production of lactate, as denoted by the sole reaction (15) (i.e., not coupled with reaction ( 14)), and excess of glutamine ends in the production of ammonia and lactate but also of other amino acids such as alanine among others as described by reactions (11) and (13).Lastly, the reaction (12) translates the consumption of glucose and glutamine for cell maintenance.The latter are two energy sources for the production of ATP and reduced pyridine nucleotides, essential for cell life.In that respect, all the chemical reactions deduced from the reduced matrix of EFMs have a biological meaning, as expected, consolidating the usefulness of the proposed reduction procedure.
The overflow metabolism phenomenon can be observed in the values of Φ i for i ∈ [1, 5] in Figure 11.Indeed, Φ 5 exhibits large values until about 100 h, before a significant decrease.This highlights a potential overflow metabolism on glucose during that first period of time, which vanishes at 100 h when glucose is almost completely depleted.This can be seen in the glucose concentration-time profile in Figure 12.A similar decrease does not appear at the levels of Φ 1 and Φ 3 .This can be interpreted as an overflow metabolism on glutamine that takes place during the whole duration of the culture as the glutamine concentration remains above a value that must be greater than the critical level (see Figure 12).

Conclusions
This study establishes a complete procedure to derive a small macroscopic bioreaction model from metabolic networks using the concept of elementary flux vectors to infer a reaction scheme and neural networks to describe the reaction rates in terms of the culture species.
When the size of the metabolic network increases, depending on the level of description required, the number of elementary modes explodes and it is important to use systematic procedures to reduce and select the most informative modes in view of establishing the corresponding minimal bioreaction scheme.For this purpose, the modular reduction procedure includes several steps: (i) the initial generation of the modes (the whole set or only a representative subset), (ii) the biological interpretation of the reaction scheme, (iii) a reduction up to Ω modes using collinearity between vectors, followed by a random elimination and (iv) the selection of a minimal set of Λ EFMs (Λ < m e ) leading to models with fewer reactions than components.
Once such a minimal model has been derived, artificial neural networks can be exploited to model the specific rates of each chemical reaction.In this study, a multilayer perceptron with one hidden layer and a small number of neurons is successfully applied to this modeling task.To further reduce the number of parameters, two simple pruning algorithms are also applied with success.Pruning algorithms appear particularly useful because of the fully connected architecture of MLPs.Further investigation of a suitable compromise between the complexity of the neural network architecture and the use of more efficient pruning algorithms could be interesting research directions.
The hybrid physical-neural models show good predictive capability in a case study related to perfusion cultures of hybridoma cells.The hybrid model is biologically interpretable in terms of reaction schemes and the reproduction of important mechanisms such as overflow metabolism.
The main advantage of the proposed procedure is the rapid development of the hybrid dynamic models (once the experimental data have been collected) whose stoichiometry and kinetics can be obtained in a systematic way.The resulting dynamic models could be exploited for process monitoring (i.e., the development of software sensors) and modelbased control, such as model predictive control.Further case studies are required to explore and consolidate these perspectives.

Figure 1 .
Figure 1.Hybrid modeling using an elementary flux mode reduction procedure and neural networks.N is the stoichiometric matrix, N e is the stoichiometric matrix of extracellular measurements, ν m is the vector of experimental uptake and excretion rates, K describes the stoichiometry of the bioreactions,

Figure 3 .
Figure 3.The modular EFM selection procedure.N is the stoichiometric matrix, N e is the stoichiometric matrix of extracellular measurements, ν m is the vector of experimental uptake and excretion rates, E is the matrix of EFMs and K describes the stoichiometry of the bioreactions.Ω is a first target number for the reduced set of EFMs, close to the number of measured extracellular species m e , and Λ is a second target number for the reduced set of EFMs, smaller than m e .NN Kinetic Modeling A perceptron, as illustrated in Figure 4, involves a set of input signals x i weighted with synaptic coefficients w i plus a bias b:

Figure 4 .
Figure 4.The structure of a perceptron: x i represent the input signals, w i denote synaptic coefficients, b is a bias, f is an activation function and y is the output signal.

Figure 5 .
Figure 5.A multilayer perceptron with one hidden layer with 3 neurons; 6 input signals in the input layer, which are the concentrations of culture species; and 5 output signals in the output layer, which are the reaction rates of each bioreaction.The biases are denoted by b i and activation functions in the hidden layer are sigmoids while activation functions in the output layer are identity functions.

Figure 6 .
Figure 6.Architecture of an MLP before and after pruning-synapses and neurons pruning.

5 Figure 9 .Figure 10 .
Figure 9.Time evolution of the specific reaction rates (in mM•h −1 ) in dataset # 2-target values (black curves); prediction with an MLP with ten neurons (red curves), five neurons (magenta curves) and three neurons (blue curves) in the hidden layer.

5 Figure 11 .
Figure11.Time evolution of the specific reaction rates (in mM•h −1 ) in dataset # 2-target values (black curves); prediction with MLP with three neurons in the hidden layer (red curves), and with three neurons and nonlinear re-identification of weights and biases (magenta curves); pruned MLP with three neurons (blue curves).

Figure 12 .
Figure12.Prediction of the time evolution of the extracellular measured species in dataset # 2numerical results from the reduction procedure before kinetic identification (black curves); MLP with three neurons in the hidden layer and nonlinear re-identification of the parameters (red curves); pruned MLP with three neurons (blue curves).

Figure 13 .
Figure13.Prediction of the time evolution of the extracellular measured species in dataset # 2pruned MLP with three neurons using magnitude pruning algorithm (blue curves); pruned MLP with three neurons using a CV pruning strategy (red curves).

Table 1 .
Summary of culture conditions for experiments in perfusion.

Table 2 .
The stoichiometric coefficients of the reaction network.

Table 3 .
Number of network parameters for six input signals and five output signals.

Table 4 .
Coefficients of variation (in %) of the estimated weights and biases-colored boxes are removed synaptic connections or biases.

Table 5 .
Results of the magnitude pruning algorithm-colored boxes are removed synaptic connections or biases.