BACKGROUND
knowledge of intracellular metabolic fluxes is crucial to understand how different pathways interact and their relative importance within the overall metabolic processes.
metabolic flux analysis  is an established methodology that allows the quantification of such intracellular fluxes.
in mfa, intracellular fluxes are calculated by applying steady-state material balances around intracellular metabolites.
in general the number of reactions exceeds the number of metabolites resulting in undetermined systems of algebraic equations  <cit> .
such systems can be solved after measurement of the missing fluxes, which are typically uptake rates of substrates and secretion rates of metabolites, and also intracellular fluxes when the former are not enough.
the determination of metabolic flux distribution in undetermined systems may also be obtained by flux-balance analysis   <cit> .
in fba, unknown fluxes are determined by linear programming , whereby a given objective function, related to a given cellular physiological state, is optimised.
typically, the maximisation of the growth flux defined in terms of biosynthetic requirements is used as the objective function  <cit> .
in general, fba provides flux distribution for a desirable physiological state, it is however uncertain that the provided solution is unique  <cit> .
frequently multiple optima are obtained which are a consequence of the existence of redundant pathways in the metabolic network conferring structural robustness to cells  <cit> .
metabolic pathway analysis  is another flux-based analysis method.
mpa, unlike fba, do not look only at the properties of solutions selected by the statement of an objective, but study the full range of achievable biochemical network states that are provided by the solution space.
network-based mpa has focused on two approaches, elementary modes  and extreme pathways   <cit> .
these approaches are very similar being eps a subset of ems.
in certain network topologies the sets of ems and eps coincide.
they are both unique for a given network and can be considered as nondecomposable steady state flux distributions using a minimal set of reactions.
the difference is that ep analysis decouples all internal reversible reactions into two separate irreversible reactions  and ems analysis accounts for reaction directionality.
in this work we have adopted the ems concept since it has broad application; eps analysis can exclude important routes of the network giving misleading results  <cit> .
mfa focuses on single flux distributions, but in a complex metabolic network of reactions there is a space of admissible flux distributions.
the mpa allows the transition from a reaction based perspective to a pathway-oriented view of metabolism because each feasible steady state flux distribution can be represented as a nonnegative combination of eps or ems  <cit> .
although flux-based analysis methods have been mainly used for metabolic engineering  <cit> , they may also be useful in other phases of the bioprocess development cycle, namely for advanced bioreactor monitoring and control  <cit> .
the ems concept is particularly attractive since it reduces network complexity to a minimal set of reactions.
provost and bastin  <cit>  exemplified the use of the ems concept for dynamic modeling of a cho culture.
the main objective of the present study is to derive an optimal control method that incorporates the knowledge of the metabolic network of the biological system under study using the ems technique.
model-based off- and on-line control techniques are today well established in both theoretical and practical terms, and have been widely applied for bioprocess optimisation .
the success of such methods is critically dependent on the quality of the supporting mathematical model.
not only accuracy in describing previously measured data but mainly the capacity to predict behaviour in unexplored states is the key for success.
in previous studies  <cit> , an iterative batch-to-batch optimisation scheme was developed and applied to the optimisation of recombinant bhk- <dig> cultures.
the method is based on the premise that in general the biological system under consideration is only partially known or even poorly known in a mechanistic sense.
following this principle, a flexible hybrid parametric/nonparametric representation of the biological system was adopted to support the batch-to-batch optimisation scheme.
it was verified that the model generalization capacity increases as more reliable mechanistic knowledge of cells is incorporated in the hybrid model.
the algorithm presented here is essentially an extension to the previous batch-to-batch optimisation scheme whereby the knowledge concerning the metabolic network is incorporated in the optimisation algorithm.
the methods will be exemplified through the application to a recombinant bhk- <dig> culture expressing the fusion glycoprotein igg1-il <dig>
RESULTS
proposed methodology
the proposed methodology for bioprocess monitoring and control is represented schematically in figures  <dig>   <dig>  the backbone of this methodology is the hybrid semiparametric model structure shown in figure  <dig>  the main design principle is flexible integration of knowledge concerning the metabolism, transport phenomena and empirical process data.
the method contemplates the possibility of missing parts of the metabolism  and of unknown reaction kinetics and underlying transduction mechanisms.
whenever knowledge is missing, empirical data modeling is called to fill the gaps.
the starting point is the establishment of the metabolic network of the biological system under study.
firstly, the metabolic network is analyzed using the elementary modes technique.
the overall network is decomposed into structural subunits, the ems, which are the simplest paths connecting substrates with end-products.
this structural analysis identifies all compounds  taken up and/or secreted to the abiotic phase, which essentially define the system state space vector.
the bioreactor dynamics are subsequently described by the material balance equations of each component occurring in the ems.
the ems kinetics are identified with data from exploratory experiments, using chemometric techniques.
once the model is properly validated, it can be used for on-line intracellular flux distribution monitoring and for on-line process performance optimisation .
so the next step is the on-line implementation of the newly developed hybrid model for process monitoring and/or optimisation in the sense of maximizing the process performance by manipulating the control inputs, i.e., the optimal control problem  <cit> .
the performance function includes a penalty term that accounts for the risk of model unreliability, i.e., extrapolation outside the model trust region.
the empirical parameters are re-estimated, followed by the re-optimisation of the future process time course whenever new measurements of the process state are performed.
in the lines below we describe in detail the steps involved.
elementary modes analysis
we consider a generic metabolic network with m metabolites and q reactions such as the network represented in figure  <dig>  assuming balanced growth and negligible dilution, the fundamental steady state mass balance equations on intracellular metabolites are expressed as follows:
{nv=0vk>0
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaadagabeqaauaabaqaceaaaeaaieqacqwfobgtcqwf2bgdcqgh9aqpcqwfwaamaeaacqwg2bgddawgaawcbagaem4aasgabeaakiabg6da+iabicdawaaaaiaawuhaaiaaxmaacawljawaaewaaeaacqaixaqmaiaawicacaglpaaaaaa@3b17@
with n = {nij} a m × q stoichiometric matrix and v = {vj} the vector of q metabolic fluxes with vj denoting the net specific rate of reaction j. some of the q reactions are irreversible, thus the respective fluxes must be nonnegative: vk >  <dig> with k denoting the irreversible reactions in the metabolic network.
the universe of solutions of system  forms a convex polyhedral cone in the solution space  <cit> .
it is a property of this system that the infinite set of solutions forming the convex polyhedral cone may be expressed as nonnegative linear combinations of a finite set of n fundamental vectors ei called elementary modes :
v=∑i=1nλiei∀v∈ℜq
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaafaqabeqacaaabagaemodaynaeyypa0zaaabcaeaaiigacqwf7oabdawgaawcbagaemyaakgabeaakiabdwgalnaabaaaleaacqwgpbqaaeqaaaqaaiabdmgapjabg2da9iabigdaxaqaaiabd6gaubqdcqghris5aagcbagaeyiaiiiaemodaynaeyici4saeyihhi8aawbaasqabeaacqwgxbqcaaaaaogaaczcaiaaxmaadaqadaqaaiabikdayagaayjkaiaawmcaaaaa@4706@
these elementary modes e = {ei} obey to constraints  and additionally to the elementarily constraint stating that there is no other non-null flux vector involving proper subsets of the reactions participating in that particular em  <cit> .
in the context of the present work, there are two most important features of the ems analysis: it allows to identify all possible pathways for the conversion of substrates into products, and opens the way to the quantification of the relative importance of pathways at a given process stage.
the non-null elements in each elementary mode, ei, define a subset of active reactions of the overall metabolic network n, which are essentially pathways for the conversion of substrates into products.
hybrid dynamical model formulation
the knowledge acquired from the ems analysis is integrated in the hybrid model structure represented in figure  <dig>  this structure allows the introduction of the a priori knowledge concerning the metabolic reactions and intracellular kinetics, but it is also open to the possibility of existing missing parts in both cases.
the ems analysis provides a stoichiometric matrix, kem, that may be augmented  if important parts of the metabolism are missing.
for instance, the energetic and biosynthetic metabolism may be known, but product metabolism may be unknown.
in this case the product concentration or other compounds and respective reactions kinetics are included through matrix i independently of the ems.
even though information about cellular components is growing rapidly, enzymes concentrations and intracellular kinetic data are difficult to obtain.
in this work the ems kinetics are defined by the product of two functions of the system state: a mechanistic/empirical function  and a nonparametric function that represents what is extracted from data without mechanistic interpretation.
this model structure can be formulated mathematically by the following two equations, which may be regarded as a general hybrid model for ideal bioreactor systems  <cit> :
dcdt=r−dc+u
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaadawcaaqaaiabdsgakhqabiab=ngajbqaaiabdsgakjabdsha0baacqgh9aqpcqwfybgccqggoaakcqwfjbwycqggsaalcqwf3bwdcqggpaqkcqghsislcqwgebarcqwfjbwycqghrawkcqwf1bqdcawljagaaczcamaabmaabagaeg4mamjaeeyyaegacagloagaayzkaaaaaa@4487@
r = k⟨ϕj × ρj⟩j =  <dig> ...,m
with c a vector of n concentrations in the liquid phase, r a vector of n volumetric reaction rates, k a n × m coefficients matrix obtained from the elementary modes analysis, ϕ are m kinetic functions established from mechanistic knowledge, ρj are m unknown kinetic functions, w a vector of parameters that must be estimated from data, d is the dilution rate, and u is a vector of n volumetric input rates .
identification of ems kinetics
the reaction rate of elementary mode j is defined by the product ϕj × ρj.
the function ϕj represents "known" mechanisms whereas ρj represents unknown mechanisms.
redundancy and degeneracy are common problems in the determination of fluxes of biological networks  <cit> .
it is very important to define a priori the conditions under which metabolic fluxes are identifiable.
a rank of matrix k equal to the number of unknown ems and an equivalent number of measured states are necessary conditions for the identifiability of system .
if identifiable, the unknown functions ρj can be extracted from data using chemometric techniques such as multilinear regression, partial least squares, artificial neural networks and many other.
in the frame of hybrid modeling, neural networks have been the most widely used technique for reaction kinetics modeling in biosystems  <cit> .
we used a backpropagation neural network with a single hidden layer to define ρj:
ρ = ρmaxs + b2)
s=11+e−x
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqqgzbwccqggoaakcqqg4baecqggpaqkcqgh9aqpdawcaaqaaiabigdaxaqaaiabigdaxiabgucariabbwgalnaacaaaleqabagaeyoei0iaeeieaghaaaaakiaaxmaacawljawaaewaaeaacqai0aancqqgibgyaiaawicacaglpaaaaaa@3e14@
with ρ = ⟨ρj⟩ a vector of m unknown reaction rates, ρmax a vector of scaling factors with dim = m, w <dig>  b <dig>  w <dig> and b <dig> are parameter matrices associated with connections between the nodes of the network, w is a vectored form of w <dig>  b <dig>  w <dig>  b <dig>  and s the sigmoid activation function.
a batch neural network training method was adopted, whereby the parameters w are estimated in the sense of least squares employing a quasi-newton optimiser with gradients calculated by the sensitivities method  <cit>  as described in the methods section.
dynamic optimisation of culture operation: optimal control problem
in the dynamic optimisation step the process performance is optimised with respect to control inputs.
this problem may be formulated mathematically as follows:
max⁡tb,c,u{j=f)+∫0tbg,u)dτ}
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaadawfqaqaaigbc2gatjabcggahjabciha4bwcbagaeeidaq3aasbaawqaaiabbkgaibqabawccqggsaalieqacqwfjbwycqggoaakcqaiwaamcqggpaqkcqggsaalcqwf1bqdcqggoaakcqqg0badcqggpaqkaeqaaowaaiwabeaacqqgkbgscqgh9aqpcqqgmbgzcqggoaakcqqgjbwycqggoaakcqqg0baddawgaawcbagaeeoyaigabeaakiabcmcapiabcmcapiabgucarmaapedabagaee4zacmaeiikagiae83yammaeiikagiaeqixdqnaeiykakiaeiilawiae8xdaunaeiikagiaeqixdqnaeiykakiaeiykakiaeeizaqmaeqixdqhaleaacqaiwaamaeaacqqg0baddawgaaadbagaeeoyaigabeaaa0gaey4kiipaaogaay5eaiaaw2haaiaaxmaacawljawaaewaaeaacqai1aqnaiaawicacaglpaaaaaa@66e5@
with j the performance index, tb the batch time, f a terminal performance function and g a time-dependent performance function.
the algorithm used was the micro-genetic algorithm  <cit>  coded by carroll  <cit> .
for simplicity, a piecewise constant approximation of the control inputs u was adopted.
the optimisation  is subject to the constraint defined by the hybrid dynamical model – ) but possibly also by other equality and inequality constraints regarding process states, cellular states and control inputs.
due to the use of nonparametric functions, namely of the neural network function , it is important to evaluate the unreliability risk during the optimisation.
after the em identification step, the measured input space is clustered by ellipsoidal functions .
the clustered input space forms the trust region, wherein the function  is considered to be reliable.
optimisation  is then further constrained by the risk of function  inputs being outside the trust region.
the technique is described in detail in the methods section.
case study: optimisation of recombinant bhk cultures
process description
to illustrate the methodology described above it will be applied to a recombinant baby hamster kidney  culture expressing a fusion glycoprotein  intended for cancer therapy  <cit> .
the experiments were carried out in serum free and protein free medium .
the batch cultures were set up in a  <dig>  <dig> reactor volume and the fed-batch cultures were set up at  <dig> different volume scales .
sparger aeration was employed.
dissolved oxygen concentration was set at 15% of air saturation.
agitation rate used was  <dig> rpm; ph was set as  <dig>  and controlled through the addition of co <dig>  experimental data of viable cells concentration and six extracellular species  were collected.
analytical techniques are described elsewhere  <cit> .
elementary modes
bhk-21a cells use glucose and glutamine as major sources of carbon and energy, and produce lactate and ammonia as toxic by-products.
a reduction of this waste production will improve both cellular growth and glycoprotein  synthesis.
figure  <dig> shows the metabolic network adopted in this work  <cit> .
as catabolic routes, the network includes the glucose and glutamine fluxes through glycolysis, glutaminolysis and tca cycle.
the amino acids metabolism was not considered; it was assumed that all of them are provided by the culture medium.
the elementary modes of the  <dig> reactions that compose the catabolism were calculated using the fluxanalyser program  <cit> .
this system has five ems, each one consisting of collections of reactions steps .
the hypothesis of balanced growth allows the elimination of the intermediate metabolites resulting in a simplified set of reactions  connecting extracellular substrates  with end-products .
the first elementary mode corresponds to the glucose flux converted into lactate; the second is the complete oxidation of glucose via tca cycle ; the third mode is the complete oxidation of glutamine  and the fourth and fifth modes are partial oxidations of glutamine in alanine and lactate, respectively.
biomass and product synthesis
in addition to the catabolism, the anabolism and product synthesis must also be considered.
for the sake of simplicity, the anabolism  was represented as a lumped equation combining the precursors of the main cellular building blocks .
the stoichiometry established by wei-shou hu and coworkers  <cit>  for an hybridoma cell line was adopted in this work.
<dig> glc +  <dig> gln +  <dig> ala +  <dig> gly +  <dig> val +  <dig> leu +  <dig> ile +  <dig> met +  <dig> pro +  <dig> phe +  <dig> try +  <dig> ser +  <dig> thr +  <dig> asn +  <dig> tyr +  <dig> cys +  <dig> lys +  <dig> arg +  <dig> his +  <dig> asp +  <dig> glu → biomass
the igg1-il <dig> synthesis was also represented as a lumped equation as follows:
<dig> gln +  <dig> ala +  <dig> gly +  <dig> val +  <dig> leu +  <dig> ile +  <dig> met +  <dig> pro +  <dig> phe +  <dig> try +  <dig> ser +  <dig> thr +  <dig> asn +  <dig> tyr +  <dig> cys +  <dig> lys +  <dig> arg +  <dig> his +  <dig> asp +  <dig> glu → igg1-il <dig>
this equation is based on the amino acid composition of the antibody igg <dig> plus the interleukin il <dig> .
the carbohydrate content of this fusion glycoprotein was neglected since the glucose contribution is extremely small compared to overall glucose consumption.
hybrid model structure
the ems analysis provides a simplified reaction mechanism based on which the following hybrid model structure  is formulated:
ddt=−d+
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaadawcaaqaaiabbsgakbqaaiabbsgakjabbsha0baadawadaqaauaabeqaheaaaaqaaiabbifaynaabaaaleaacqqg2bgdaeqaaagcbagaee4rackaeeibawmaee4yamgabagaee4rackaeeibawmaeeoba4gabagaeeitawkaeeyyaemaee4yamgabagaeeyqaekaeeyba0maeeyba0gabagaeeyqaekaeeibawmaeeyyaegabagaeeysakkaee4zacmaee4raceaaagaay5waiaaw2faaiabg2da9maadmaabaqbaeqabcgbaaaaaeaacqaiwaamaeaacqaiwaamaeaacqaiwaamaeaacqaiwaamaeaacqaixaqmaeaacqaiwaamaeaacqghsislcqaixaqmaeaacqghsislcqaixaqmaeaacqaiwaamaeaacqaiwaamaeaacqghsislcqaiwaamcqgguaglcqaiwaamcqaiyagmcqaiwaamcqai4aaoaeaacqaiwaamaeaacqaiwaamaeaacqaiwaamaeaacqghsislcqaixaqmaeaacqghsislcqaixaqmaeaacqghsislcqaiwaamcqgguaglcqaiwaamcqaizawmcqai3awncqai3awnaeaacqghsislcqaiwaamcqgguaglcqaiwaamcqaixaqmcqaiwaamcqai0aanaeaacqaiyagmaeaacqaiwaamaeaacqaiwaamaeaacqaiwaamaeaacqaiwaamaeaacqaiwaamaeaacqaiwaamaeaacqaiwaamaeaacqaiyagmaeaacqaixaqmaeaacqaiwaamaeaacqaiwaamaeaacqaiwaamaeaacqaiwaamaeaacqaiwaamaeaacqaixaqmaeaacqghsislcqaiwaamcqgguaglcqaiwaamcqaixaqmcqaizawmcqaizawmaeaacqghsislcqaiwaamcqgguaglcqaiwaamcqaixaqmcqaixaqmcqaiyagmaeaacqaiwaamaeaacqaiwaamaeaacqaiwaamaeaacqaiwaamaeaacqaiwaamaeaacqaixaqmaaaacaglbbgaayzxaawaamwaaeaafaqabegbbaaaaeaacqahbpgcdawgaawcbagaegymaedabeaakiabcicaoiabbeeahjabbygasjabbogajjabbycasiabbccagiabbeeahjabbygasjabb6gaujabbycasiabbccagiabbgeabjabb2gatjabb2gatjabcmcapiabgwsixlabcicaoiabbifaynaabaaaleaacqqgwbgvaeqaaogaee4rackaeeibawmaee4yammaeiykakcabagaeqywdi3aasbaasqaaiabikdayaqabagccqggoaakcqqghbwrcqqgsbabcqqgjbwycqqgsaalcqqggaaicqqghbwrcqqgsbabcqqgubgbcqqgsaalcqqggaaicqqgbbqqcqqgtbqbcqqgtbqbcqggpaqkcqghfly1cqggoaakcqqgybawdawgaawcbagaeeovayfabeaakiabbeeahjabbygasjabbogajjabcmcapaqaaiabeg8aynaabaaaleaacqaizawmaeqaaogaeiikagiaee4rackaeeibawmaee4yammaeeilawiaeeiiaaiaee4rackaeeibawmaeeoba4maeeilawiaeeiiaaiaeeyqaekaeeyba0maeeyba0maeiykakiaeyyxictaeiikagiaeeiwag1aasbaasqaaiabbafawbqabagccqqghbwrcqqgsbabcqqgubgbcqggpaqkaeaacqahbpgcdawgaawcbagaeginaqdabeaakiabcicaoiabbeeahjabbygasjabbogajjabbycasiabbccagiabbeeahjabbygasjabb6gaujabbycasiabbccagiabbgeabjabb2gatjabb2gatjabcmcapiabgwsixlabcicaoiabbifaynaabaaaleaacqqgwbgvaeqaaogaee4rackaeeibawmaeeoba4maeiykakcabagaeqivd0maeiikagiaee4rackaeeibawmaee4yammaeeilawiaeeiiaaiaee4rackaeeibawmaeeoba4maeeilawiaeeiiaaiaeeyqaekaeeyba0maeeyba0maeiykakiaeyyxictaeiikagiaeeiwag1aasbaasqaaiabbafawbqabagccqggpaqkaeaacqahapaccqggoaakcqqghbwrcqqgsbabcqqgjbwycqqgsaalcqqggaaicqqghbwrcqqgsbabcqqgubgbcqqgsaalcqqggaaicqqgbbqqcqqgtbqbcqqgtbqbcqggpaqkcqghfly1cqggoaakcqqgybawdawgaawcbagaeeovayfabeaakiabbeeahjabbygasjabb6gaujabcmcapaaaaiaawufacagldbaacqghsislcqqgebardawadaqaauaabeqaheaaaaqaaiabbifaynaabaaaleaacqqg2bgdaeqaaagcbagaee4rackaeeibawmaee4yamgabagaee4rackaeeibawmaeeoba4gabagaeeitawkaeeyyaemaee4yamgabagaeeyqaekaeeyba0maeeyba0gabagaeeyqaekaeeibawmaeeyyaegabagaeeysakkaee4zacmaee4raceaaagaay5waiaaw2faaiabgucarmaadmaabaqbaeqabcqaaaaabagaegimaadabagaeeoray0aasbaasqaaiabbeeahjabbygasjabbogajbqabaaakeaacqqggbgrdawgaawcbagaee4rackaeeibawmaeeoba4gabeaaaoqaaiabicdawaqaaiabicdawaqaaiabicdawaqaaiabicdawaaaaiaawufacagldbaacawljagaaczcamaabmaabagaegyoakdacagloagaayzkaaaaaa@76dc@
the state space vector is formed by the concentrations of compounds that intervene in the final reactions set  and additionally the concentrations of viable cells, xv, and product, igg:
c = t.
carbon dioxide was excluded because its concentration was not measured and because it doesn't interfere with the dynamics of the remaining variables .
the stoichiometric matrix, k, is established from the elementary modes of table  <dig>  but it also accounts for cell growth  and product formation  as lumped equations of glucose, glutamine and alanine.
it should be noted that the 5th em was not included in k because preliminary results showed that this em has negligible flux.
this observation is in agreement with some published works  <cit> , stating that lactate is mainly produced from glucose, being the percentage coming from glutamine very low .
the volumetric reaction rates of the ems were defined on the basis of the following assumptions:
all reaction rates are specific ,
the metabolic reactions considered are all irreversible  and therefore the respective reaction rates are nonnegative
uncertainty in relation to kinetic constants and possible unknown saturation and inhibition effects.
only the concentrations of glucose, glutamine and ammonia have a significant impact on the specific reaction kinetics  <cit> .
lactate never reaches inhibitory levels in our experiments.
in the reaction rates of eq.
the term in parenthesis represents the a priori knowledge concerning the kinetics of the particular reaction  and ) whereas the ρi, μ and π terms represent the uncertainty concerning the reaction kinetics ) and are functions of three state variables ).
with this particular formulation, the vector of known kinetic functions is given by:
ϕ = t,
whereas the vector of unknown kinetics is given by:
ρ = t = ρ
the last term in eq.
is the control input vector u that in our case accounts for the volumetric feeding of glucose, fglc, and glutamine, fgln.
identification of the em kinetics
an important point for the identification of unknown flux functions  is that the rank of k is  <dig>  thus the measurement of   =  <dig> > rank) is sufficient for the observability of the em kinetics.
the other relevant point is the availability of sufficiently "rich" measurements to identify the "true" fluxes.
preliminary simulation tests showed that, for the system structure of eq , the "true" fluxes can be identified under typical fed-batch conditions .
experimental data of seven experiments  were used for the identification of the em kinetics.
data of  <dig> experiments were used for parameter calibration whereas data of  <dig> experiments were used for model validation.
the concentrations in the state space vector  were analyzed off-line according to methods published elsewhere  <cit> .
the neural network has three inputs, t and  <dig> outputs as defined by eq.
.
the number of hidden nodes was tuned heuristically in the sense of minimizing the error of the validation data set.
the best result was obtained with  <dig> hidden nodes giving a total number of network parameters equal to  <dig>  the output scaling factors reflect the maximum kinetic rates and were defined as ρmax = t.
the hybrid modeling results in terms of measured and predicted state variables are presented in figure  <dig> for both training and validation data sets.
examples of em kinetics identification are provided in figure  <dig>  the hybrid model was able to describe simultaneously all seven experiments with acceptable accuracy.
in particular, the results with the validation data set strengthen the predictive potential of the model.
metabolic interpretation
in general, the metabolic activity of bhk cells during the cell growth phase is higher in the batch culture than in the fed-batch culture.
the em fluxes in the fed-batch experiments seem to be much more controlled than in the batch experiments.
the high levels of glucose and glutamine at the beginning of the batch culture are mostly directed toward the overflow metabolism, i.e., waste production of lactate and alanine .
the fed-batch culture, which had lower glucose and glutamine concentrations during the cell growth phase, started with substantially lower consumption rates of these nutrients, while maintaining the flux to biomass.
glucose is consumed for biomass synthesis and is metabolized through elementary modes em <dig> and em <dig>  the consumption of this nutrient differs significantly from batch to fed-batch cultures.
cells use glucose in a more efficient way in fed-batch than in batch cultures particularly at the beginning of the culture, since the glucose metabolized to lactate  is much higher in the batch experiment.
these results are in agreement with several published works for other mammalian cells  <cit> , where it is reported that high glycolytic activity of animal cells results from high residual glucose.
glutamine is consumed for product and biomass synthesis and is metabolized through elementary modes em <dig> and em <dig>  the most energetic pathway involving glutamine, em <dig>  is practically inactive during the cell growth phase of the batch culture, being glutamine preferentially converted into alanine.
on the opposite, this elementary flux mode is an important pathway in the course of the fed-batch culture, representing 50% of the total glutamine consumption.
the higher production rate of alanine at high glutamine levels  is consistent with observations made by doverskog et al.
<cit>  and vriezen et al.
<cit> .
concerning the specific formation rate of glycoprotein, πigg, the product formation is consistently more stable in the fed-batch culture than in the batch culture.
the product synthesis rate oscillates between  <dig>  and  <dig>  mg 10- <dig> cells h- <dig> in the former case whereas in the latter case the reaction is much lower in the beginning.
this appears to be correlated with the overflow metabolism in the batch experiments, which seems to be detrimental for product synthesis.
on-line culture optimisation
the hybrid model was further used for on-line optimisation of a fed-batch bhk culture.
the model parameters were re-tuned on-line using the data from the running cultivation.
a batch training scheme was adopted, whereby the data of the running cultivation along with the data of historical experiments were used for model adjustment .
some variables, namely glucose, glutamine, lactate and viable cells concentration, can be measured off-line, with the results available after a short period of time .
therefore, at sampling times of  <dig> to  <dig> h, these measurements were stored in the training data set and then used for parameter identification using the same strategy previously described in the identification of the em kinetics section.
the only difference was that the initial parameter values were those obtained in the off-line training procedure.
the ammonia, alanine and igg concentrations were exclude form this model adjustment step because they were quantified only at the end of the experiment.
after the parameters retuning step, the glucose and glutamine feeding rates were re-optimised  in the sense of maximising the total amount of igg1-il <dig> produced at the end of the experiment .
max⁡uj=cigg−il2vequation 0≤fglu≤fglu,max0≤fgln≤fgln,max−1∫t0tfriskdτ≤riskmax⁡
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaafaqaaeqbcaaaaeaadawfqaqaaigbc2gatjabcggahjabciha4bwcbagaemydauhabeaaaoqaaiabdqeakjabg2da9iabdoeadnaabaaaleaacqwgjbqscqwgnbwzcqwghbwrcqghsislcqwgjbqscqwgmbatcqaiyagmaeqaaogaeiikagiaemidaq3aasbaasqaaiabdagambqabagccqggpaqkcqwgwbgvcqggoaakcqwg0baddawgaawcbagaemozaygabeaakiabcmcapaqaaaqaaiabbwgaljabbghaxjabbwha1jabbggahjabbsha0jabbmgapjabb+gavjabb6gaujabbccagmaabmaabagaegyoakdacagloagaayzkaaaabaaabagaegimaajaeyizimqaemoray0aasbaasqaaiabdeganjabdygasjabdwha1bqabagccqghkjyocqwggbgrdawgaawcbagaem4zacmaemibawmaemydaunaeiilawccbigae8xba0mae8xyaemae8heaghabeaaaoqaaaqaaiabicdawiabgsmijkabdaeagnaabaaaleaacqwgnbwzcqwfsbabcqwfubgbaeqaaogaeyizimqaemoray0aasbaasqaaiabdeganjab=xgasjab=5gaujabcycasiab=1gatjab=fgahjab=hha4bqabaaakeaaaeaacqggoaakcqwg0baddawgaawcbagaemozaygabeaakiabgkhitiabdsha0naabaaaleaacqaiwaamaeqaaogaeiykakyaawbaasqabeaacqghsislcqaixaqmaagcdawdxbqaaiabdkfasjabdmeajjabdofatjabduealnaadmaabaacbegae43yammaeiikagcccigae0hxdqnaeiykakcacaglbbgaayzxaagaemizaqmae0hxdqnaeyizimqaemouailaemysakkaem4uamlaem4sas0aasbaasqaaigbc2gatjabcggahjabciha4bqabaaabagaemidaq3aasbaawqaaiabicdawaqabaaaleaacqwg0baddawgaaadbagaemozaygabeaaa0gaey4kiipaaagccawljagaaczcamaabmaabagaegymaejaeg4mamdacagloagaayzkaaaaaa@aed9@
the optimisation  is subjected to the constraints of the process dynamics , upper and lower bounds for the glucose and glutamine feeding rates, and the maximum risk of model unreliability, riskmax .
the increase in volume is negligible in this problem, thus it was not considered in optimisation .
the risk constraint states that the average risk must not exceed a given maximum level defined by the user.
this restricts the feasibility domain to low risk regions and is essential for process optimisation supported by hybrid models since the black-box model  predictions degrade in regions of the input space with sparse measurements.
the micro-genetic algorithm coded by  <cit>  was used to solve optimisation .
the population size and number of generations was  <dig> and  <dig>  respectively.
a maximum risk level of 35%  was adopted during the on-line optimisation experiment.
at each sampling time, the flow rates of the feeding pumps were updated according to the re-optimised feeding profiles of glucose and glutamine.
medium optimisation
higher productivities are likely if the initial medium composition is optimised.
the optimisation results of initial glucose and glutamine concentrations along with the corresponding feeding strategies are shown in figure 8a.
the medium should have initially low levels of glucose and glutamine.
as for the feeding strategies, low levels of glucose and glutamine are optimal during the cell growth phase , whereas during the cell death phase, glutamine should slightly increase and glucose should decrease to a concentration of  <dig>  mm.
a significant increase of product titre was predicted under such cultivation conditions.
the kinetics of the elementary modes for the optimal strategy is presented in figure 8b.
the optimal strategy can be interpreted by means of intracellular flux distribution.
the results were not much different from the elementary flux distribution of the fed-batch culture .
the main difference arises from the fact that the specific product formation  increases during the cell death phase as a consequence of the increase of the glutamine concentration.
CONCLUSIONS
the current work presents a hybrid modeling approach for bioprocesses that integrates information concerning intracellular metabolic fluxes, bioreactor transport phenomena and measured data.
first, elementary flux analysis is used to reduce the biosystem metabolic network into a set of macroscopic reactions relating extracellular components only.
a bioreactor dynamical model is then established from material balance equations of the compounds which intervene in the final reaction set.
nonparametric techniques are used for identification of the elementary modes kinetics from measured data.
the method was successfully applied to a recombinant bhk-21a cell line producing a fusion glycoprotein.
a significant result is the achievement of the flux distribution over the runtime of a bioprocess.
the so obtained information can be used to identify conditions that favour product formation.
a fed-batch bhk culture performed with on-line optimisation supported by the proposed methodology allowed a 10% increase in the final productivity.
higher productivities are expected if starting nutrients concentrations are optimised.
the developed tool promises to be advantageous for optimising the productivity of fed-batch biochemical processes since the transfer and adaptation to different cell lines is reasonably straightforward

 knowledge of intracellular metabolic fluxes is crucial to understand how different pathways interact and their relative importance within the overall metabolic processes. metabolic flux analysis is an established methodology that allows the quantification of such intracellular fluxes. in mfa, intracellular fluxes are calculated by applying steady-state material balances around intracellular metabolites. in general the number of reactions exceeds the number of metabolites resulting in undetermined systems of algebraic equations. such systems can be solved after measurement of the missing fluxes, which are typically uptake rates of substrates and secretion rates of metabolites, and also intracellular fluxes when the former are not enough. the determination of metabolic flux distribution in undetermined systems may also be obtained by flux-balance analysis. in fba, unknown fluxes are determined by linear programming, whereby a given objective function, related to a given cellular physiological state, is optimised. typically, the maximisation of the growth flux defined in terms of biosynthetic requirements is used as the objective function. in general, fba provides flux distribution for a desirable physiological state, it is however uncertain that the provided solution is unique. frequently multiple optima are obtained which are a consequence of the existence of redundant pathways in the metabolic network conferring structural robustness to cells. metabolic pathway analysis is another flux-based analysis method. mpa, unlike fba, do not look only at the properties of solutions selected by the statement of an objective, but study the full range of achievable biochemical network states that are provided by the solution space. network-based mpa has focused on two approaches, elementary modes and extreme pathways. these approaches are very similar being eps a subset of ems. in certain network topologies the sets of ems and eps coincide. they are both unique for a given network and can be considered as nondecomposable steady state flux distributions using a minimal set of reactions. the difference is that ep analysis decouples all internal reversible reactions into two separate irreversible reactions and ems analysis accounts for reaction directionality. in this work we have adopted the ems concept since it has broad application; eps analysis can exclude important routes of the network giving misleading results. mfa focuses on single flux distributions, but in a complex metabolic network of reactions there is a space of admissible flux distributions. the mpa allows the transition from a reaction based perspective to a pathway-oriented view of metabolism because each feasible steady state flux distribution can be represented as a nonnegative combination of eps or ems. although flux-based analysis methods have been mainly used for metabolic engineering, they may also be useful in other phases of the bioprocess development cycle, namely for advanced bioreactor monitoring and control. the ems concept is particularly attractive since it reduces network complexity to a minimal set of reactions. provost and bastin exemplified the use of the ems concept for dynamic modeling of a cho culture. the main objective of the present study is to derive an optimal control method that incorporates the knowledge of the metabolic network of the biological system under study using the ems technique. model-based off- and on-line control techniques are today well established in both theoretical and practical terms, and have been widely applied for bioprocess optimisation. the success of such methods is critically dependent on the quality of the supporting mathematical model. not only accuracy in describing previously measured data but mainly the capacity to predict behaviour in unexplored states is the key for success. in previous studies, an iterative batch-to-batch optimisation scheme was developed and applied to the optimisation of recombinant bhk- cultures. the method is based on the premise that in general the biological system under consideration is only partially known or even poorly known in a mechanistic sense. following this principle, a flexible hybrid parametric/nonparametric representation of the biological system was adopted to support the batch-to-batch optimisation scheme. it was verified that the model generalization capacity increases as more reliable mechanistic knowledge of cells is incorporated in the hybrid model. the algorithm presented here is essentially an extension to the previous batch-to-batch optimisation scheme whereby the knowledge concerning the metabolic network is incorporated in the optimisation algorithm. the methods will be exemplified through the application to a recombinant bhk- culture expressing the fusion glycoprotein igg1-il  proposed methodology the proposed methodology for bioprocess monitoring and control is represented schematically in figures the backbone of this methodology is the hybrid semiparametric model structure shown in figure the main design principle is flexible integration of knowledge concerning the metabolism, transport phenomena and empirical process data. the method contemplates the possibility of missing parts of the metabolism and of unknown reaction kinetics and underlying transduction mechanisms. whenever knowledge is missing, empirical data modeling is called to fill the gaps. the starting point is the establishment of the metabolic network of the biological system under study. firstly, the metabolic network is analyzed using the elementary modes technique. the overall network is decomposed into structural subunits, the ems, which are the simplest paths connecting substrates with end-products. this structural analysis identifies all compounds taken up and/or secreted to the abiotic phase, which essentially define the system state space vector. the bioreactor dynamics are subsequently described by the material balance equations of each component occurring in the ems. the ems kinetics are identified with data from exploratory experiments, using chemometric techniques. once the model is properly validated, it can be used for on-line intracellular flux distribution monitoring and for on-line process performance optimisation. so the next step is the on-line implementation of the newly developed hybrid model for process monitoring and/or optimisation in the sense of maximizing the process performance by manipulating the control inputs, i.e., the optimal control problem. the performance function includes a penalty term that accounts for the risk of model unreliability, i.e., extrapolation outside the model trust region. the empirical parameters are re-estimated, followed by the re-optimisation of the future process time course whenever new measurements of the process state are performed. in the lines below we describe in detail the steps involved. elementary modes analysis we consider a generic metabolic network with m metabolites and q reactions such as the network represented in figure assuming balanced growth and negligible dilution, the fundamental steady state mass balance equations on intracellular metabolites are expressed as follows: {nv=0vk>0      mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaadagabeqaauaabaqaceaaaeaaieqacqwfobgtcqwf2bgdcqgh9aqpcqwfwaamaeaacqwg2bgddawgaawcbagaem4aasgabeaakiabg6da+iabicdawaaaaiaawuhaaiaaxmaacawljawaaewaaeaacqaixaqmaiaawicacaglpaaaaaa@3b17@ with n = {nij} a m × q stoichiometric matrix and v = {vj} the vector of q metabolic fluxes with vj denoting the net specific rate of reaction j. some of the q reactions are irreversible, thus the respective fluxes must be nonnegative: vk > with k denoting the irreversible reactions in the metabolic network. the universe of solutions of system forms a convex polyhedral cone in the solution space. it is a property of this system that the infinite set of solutions forming the convex polyhedral cone may be expressed as nonnegative linear combinations of a finite set of n fundamental vectors ei called elementary modes: v=∑i=1nλiei∀v∈ℜq      mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaafaqabeqacaaabagaemodaynaeyypa0zaaabcaeaaiigacqwf7oabdawgaawcbagaemyaakgabeaakiabdwgalnaabaaaleaacqwgpbqaaeqaaaqaaiabdmgapjabg2da9iabigdaxaqaaiabd6gaubqdcqghris5aagcbagaeyiaiiiaemodaynaeyici4saeyihhi8aawbaasqabeaacqwgxbqcaaaaaogaaczcaiaaxmaadaqadaqaaiabikdayagaayjkaiaawmcaaaaa@4706@ these elementary modes e = {ei} obey to constraints and additionally to the elementarily constraint stating that there is no other non-null flux vector involving proper subsets of the reactions participating in that particular em. in the context of the present work, there are two most important features of the ems analysis: it allows to identify all possible pathways for the conversion of substrates into products, and opens the way to the quantification of the relative importance of pathways at a given process stage. the non-null elements in each elementary mode, ei, define a subset of active reactions of the overall metabolic network n, which are essentially pathways for the conversion of substrates into products. hybrid dynamical model formulation the knowledge acquired from the ems analysis is integrated in the hybrid model structure represented in figure this structure allows the introduction of the a priori knowledge concerning the metabolic reactions and intracellular kinetics, but it is also open to the possibility of existing missing parts in both cases. the ems analysis provides a stoichiometric matrix, kem, that may be augmented if important parts of the metabolism are missing. for instance, the energetic and biosynthetic metabolism may be known, but product metabolism may be unknown. in this case the product concentration or other compounds and respective reactions kinetics are included through matrix i independently of the ems. even though information about cellular components is growing rapidly, enzymes concentrations and intracellular kinetic data are difficult to obtain. in this work the ems kinetics are defined by the product of two functions of the system state: a mechanistic/empirical function and a nonparametric function that represents what is extracted from data without mechanistic interpretation. this model structure can be formulated mathematically by the following two equations, which may be regarded as a general hybrid model for ideal bioreactor systems: dcdt=r−dc+u      mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaadawcaaqaaiabdsgakhqabiab=ngajbqaaiabdsgakjabdsha0baacqgh9aqpcqwfybgccqggoaakcqwfjbwycqggsaalcqwf3bwdcqggpaqkcqghsislcqwgebarcqwfjbwycqghrawkcqwf1bqdcawljagaaczcamaabmaabagaeg4mamjaeeyyaegacagloagaayzkaaaaaa@4487@ r = k⟨ϕj × ρj⟩j =...,m     with c a vector of n concentrations in the liquid phase, r a vector of n volumetric reaction rates, k a n × m coefficients matrix obtained from the elementary modes analysis, ϕ are m kinetic functions established from mechanistic knowledge, ρj are m unknown kinetic functions, w a vector of parameters that must be estimated from data, d is the dilution rate, and u is a vector of n volumetric input rates. identification of ems kinetics the reaction rate of elementary mode j is defined by the product ϕj × ρj. the function ϕj represents "known" mechanisms whereas ρj represents unknown mechanisms. redundancy and degeneracy are common problems in the determination of fluxes of biological networks. it is very important to define a priori the conditions under which metabolic fluxes are identifiable. a rank of matrix k equal to the number of unknown ems and an equivalent number of measured states are necessary conditions for the identifiability of system. if identifiable, the unknown functions ρj can be extracted from data using chemometric techniques such as multilinear regression, partial least squares, artificial neural networks and many other. in the frame of hybrid modeling, neural networks have been the most widely used technique for reaction kinetics modeling in biosystems. we used a backpropagation neural network with a single hidden layer to define ρj: ρ = ρmaxs + b2)     s=11+e−x      mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqqgzbwccqggoaakcqqg4baecqggpaqkcqgh9aqpdawcaaqaaiabigdaxaqaaiabigdaxiabgucariabbwgalnaacaaaleqabagaeyoei0iaeeieaghaaaaakiaaxmaacawljawaaewaaeaacqai0aancqqgibgyaiaawicacaglpaaaaaa@3e14@ with ρ = ⟨ρj⟩ a vector of m unknown reaction rates, ρmax a vector of scaling factors with dim = m, w b w and b are parameter matrices associated with connections between the nodes of the network, w is a vectored form of w b w b and s the sigmoid activation function. a batch neural network training method was adopted, whereby the parameters w are estimated in the sense of least squares employing a quasi-newton optimiser with gradients calculated by the sensitivities method as described in the methods section. dynamic optimisation of culture operation: optimal control problem in the dynamic optimisation step the process performance is optimised with respect to control inputs. this problem may be formulated mathematically as follows: max⁡tb,c,u{j=f)+∫0tbg,u)dτ}      mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaadawfqaqaaigbc2gatjabcggahjabciha4bwcbagaeeidaq3aasbaawqaaiabbkgaibqabawccqggsaalieqacqwfjbwycqggoaakcqaiwaamcqggpaqkcqggsaalcqwf1bqdcqggoaakcqqg0badcqggpaqkaeqaaowaaiwabeaacqqgkbgscqgh9aqpcqqgmbgzcqggoaakcqqgjbwycqggoaakcqqg0baddawgaawcbagaeeoyaigabeaakiabcmcapiabcmcapiabgucarmaapedabagaee4zacmaeiikagiae83yammaeiikagiaeqixdqnaeiykakiaeiilawiae8xdaunaeiikagiaeqixdqnaeiykakiaeiykakiaeeizaqmaeqixdqhaleaacqaiwaamaeaacqqg0baddawgaaadbagaeeoyaigabeaaa0gaey4kiipaaogaay5eaiaaw2haaiaaxmaacawljawaaewaaeaacqai1aqnaiaawicacaglpaaaaaa@66e5@ with j the performance index, tb the batch time, f a terminal performance function and g a time-dependent performance function. the algorithm used was the micro-genetic algorithm coded by carroll. for simplicity, a piecewise constant approximation of the control inputs u was adopted. the optimisation is subject to the constraint defined by the hybrid dynamical model – ) but possibly also by other equality and inequality constraints regarding process states, cellular states and control inputs. due to the use of nonparametric functions, namely of the neural network function, it is important to evaluate the unreliability risk during the optimisation. after the em identification step, the measured input space is clustered by ellipsoidal functions. the clustered input space forms the trust region, wherein the function is considered to be reliable. optimisation is then further constrained by the risk of function inputs being outside the trust region. the technique is described in detail in the methods section. case study: optimisation of recombinant bhk cultures process description to illustrate the methodology described above it will be applied to a recombinant baby hamster kidney culture expressing a fusion glycoprotein intended for cancer therapy. the experiments were carried out in serum free and protein free medium. the batch cultures were set up in a reactor volume and the fed-batch cultures were set up at different volume scales. sparger aeration was employed. dissolved oxygen concentration was set at 15% of air saturation. agitation rate used was rpm; ph was set as and controlled through the addition of co experimental data of viable cells concentration and six extracellular species were collected. analytical techniques are described elsewhere. elementary modes bhk-21a cells use glucose and glutamine as major sources of carbon and energy, and produce lactate and ammonia as toxic by-products. a reduction of this waste production will improve both cellular growth and glycoprotein synthesis. figure shows the metabolic network adopted in this work. as catabolic routes, the network includes the glucose and glutamine fluxes through glycolysis, glutaminolysis and tca cycle. the amino acids metabolism was not considered; it was assumed that all of them are provided by the culture medium. the elementary modes of the reactions that compose the catabolism were calculated using the fluxanalyser program. this system has five ems, each one consisting of collections of reactions steps. the hypothesis of balanced growth allows the elimination of the intermediate metabolites resulting in a simplified set of reactions connecting extracellular substrates with end-products. the first elementary mode corresponds to the glucose flux converted into lactate; the second is the complete oxidation of glucose via tca cycle; the third mode is the complete oxidation of glutamine and the fourth and fifth modes are partial oxidations of glutamine in alanine and lactate, respectively. biomass and product synthesis in addition to the catabolism, the anabolism and product synthesis must also be considered. for the sake of simplicity, the anabolism was represented as a lumped equation combining the precursors of the main cellular building blocks. the stoichiometry established by wei-shou hu and coworkers for an hybridoma cell line was adopted in this work. glc + gln + ala + gly + val + leu + ile + met + pro + phe + try + ser + thr + asn + tyr + cys + lys + arg + his + asp + glu → biomass     the igg1-il synthesis was also represented as a lumped equation as follows: gln + ala + gly + val + leu + ile + met + pro + phe + try + ser + thr + asn + tyr + cys + lys + arg + his + asp + glu → igg1-il     this equation is based on the amino acid composition of the antibody igg plus the interleukin il. the carbohydrate content of this fusion glycoprotein was neglected since the glucose contribution is extremely small compared to overall glucose consumption. hybrid model structure the ems analysis provides a simplified reaction mechanism based on which the following hybrid model structure is formulated: ddt=−d+      mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaadawcaaqaaiabbsgakbqaaiabbsgakjabbsha0baadawadaqaauaabeqaheaaaaqaaiabbifaynaabaaaleaacqqg2bgdaeqaaagcbagaee4rackaeeibawmaee4yamgabagaee4rackaeeibawmaeeoba4gabagaeeitawkaeeyyaemaee4yamgabagaeeyqaekaeeyba0maeeyba0gabagaeeyqaekaeeibawmaeeyyaegabagaeeysakkaee4zacmaee4raceaaagaay5waiaaw2faaiabg2da9maadmaabaqbaeqabcgbaaaaaeaacqaiwaamaeaacqaiwaamaeaacqaiwaamaeaacqaiwaamaeaacqaixaqmaeaacqaiwaamaeaacqghsislcqaixaqmaeaacqghsislcqaixaqmaeaacqaiwaamaeaacqaiwaamaeaacqghsislcqaiwaamcqgguaglcqaiwaamcqaiyagmcqaiwaamcqai4aaoaeaacqaiwaamaeaacqaiwaamaeaacqaiwaamaeaacqghsislcqaixaqmaeaacqghsislcqaixaqmaeaacqghsislcqaiwaamcqgguaglcqaiwaamcqaizawmcqai3awncqai3awnaeaacqghsislcqaiwaamcqgguaglcqaiwaamcqaixaqmcqaiwaamcqai0aanaeaacqaiyagmaeaacqaiwaamaeaacqaiwaamaeaacqaiwaamaeaacqaiwaamaeaacqaiwaamaeaacqaiwaamaeaacqaiwaamaeaacqaiyagmaeaacqaixaqmaeaacqaiwaamaeaacqaiwaamaeaacqaiwaamaeaacqaiwaamaeaacqaiwaamaeaacqaixaqmaeaacqghsislcqaiwaamcqgguaglcqaiwaamcqaixaqmcqaizawmcqaizawmaeaacqghsislcqaiwaamcqgguaglcqaiwaamcqaixaqmcqaixaqmcqaiyagmaeaacqaiwaamaeaacqaiwaamaeaacqaiwaamaeaacqaiwaamaeaacqaiwaamaeaacqaixaqmaaaacaglbbgaayzxaawaamwaaeaafaqabegbbaaaaeaacqahbpgcdawgaawcbagaegymaedabeaakiabcicaoiabbeeahjabbygasjabbogajjabbycasiabbccagiabbeeahjabbygasjabb6gaujabbycasiabbccagiabbgeabjabb2gatjabb2gatjabcmcapiabgwsixlabcicaoiabbifaynaabaaaleaacqqgwbgvaeqaaogaee4rackaeeibawmaee4yammaeiykakcabagaeqywdi3aasbaasqaaiabikdayaqabagccqggoaakcqqghbwrcqqgsbabcqqgjbwycqqgsaalcqqggaaicqqghbwrcqqgsbabcqqgubgbcqqgsaalcqqggaaicqqgbbqqcqqgtbqbcqqgtbqbcqggpaqkcqghfly1cqggoaakcqqgybawdawgaawcbagaeeovayfabeaakiabbeeahjabbygasjabbogajjabcmcapaqaaiabeg8aynaabaaaleaacqaizawmaeqaaogaeiikagiaee4rackaeeibawmaee4yammaeeilawiaeeiiaaiaee4rackaeeibawmaeeoba4maeeilawiaeeiiaaiaeeyqaekaeeyba0maeeyba0maeiykakiaeyyxictaeiikagiaeeiwag1aasbaasqaaiabbafawbqabagccqqghbwrcqqgsbabcqqgubgbcqggpaqkaeaacqahbpgcdawgaawcbagaeginaqdabeaakiabcicaoiabbeeahjabbygasjabbogajjabbycasiabbccagiabbeeahjabbygasjabb6gaujabbycasiabbccagiabbgeabjabb2gatjabb2gatjabcmcapiabgwsixlabcicaoiabbifaynaabaaaleaacqqgwbgvaeqaaogaee4rackaeeibawmaeeoba4maeiykakcabagaeqivd0maeiikagiaee4rackaeeibawmaee4yammaeeilawiaeeiiaaiaee4rackaeeibawmaeeoba4maeeilawiaeeiiaaiaeeyqaekaeeyba0maeeyba0maeiykakiaeyyxictaeiikagiaeeiwag1aasbaasqaaiabbafawbqabagccqggpaqkaeaacqahapaccqggoaakcqqghbwrcqqgsbabcqqgjbwycqqgsaalcqqggaaicqqghbwrcqqgsbabcqqgubgbcqqgsaalcqqggaaicqqgbbqqcqqgtbqbcqqgtbqbcqggpaqkcqghfly1cqggoaakcqqgybawdawgaawcbagaeeovayfabeaakiabbeeahjabbygasjabb6gaujabcmcapaaaaiaawufacagldbaacqghsislcqqgebardawadaqaauaabeqaheaaaaqaaiabbifaynaabaaaleaacqqg2bgdaeqaaagcbagaee4rackaeeibawmaee4yamgabagaee4rackaeeibawmaeeoba4gabagaeeitawkaeeyyaemaee4yamgabagaeeyqaekaeeyba0maeeyba0gabagaeeyqaekaeeibawmaeeyyaegabagaeeysakkaee4zacmaee4raceaaagaay5waiaaw2faaiabgucarmaadmaabaqbaeqabcqaaaaabagaegimaadabagaeeoray0aasbaasqaaiabbeeahjabbygasjabbogajbqabaaakeaacqqggbgrdawgaawcbagaee4rackaeeibawmaeeoba4gabeaaaoqaaiabicdawaqaaiabicdawaqaaiabicdawaqaaiabicdawaaaaiaawufacagldbaacawljagaaczcamaabmaabagaegyoakdacagloagaayzkaaaaaa@76dc@ the state space vector is formed by the concentrations of compounds that intervene in the final reactions set and additionally the concentrations of viable cells, xv, and product, igg: c = t.     carbon dioxide was excluded because its concentration was not measured and because it doesn't interfere with the dynamics of the remaining variables. the stoichiometric matrix, k, is established from the elementary modes of table but it also accounts for cell growth and product formation as lumped equations of glucose, glutamine and alanine. it should be noted that the 5th em was not included in k because preliminary results showed that this em has negligible flux. this observation is in agreement with some published works, stating that lactate is mainly produced from glucose, being the percentage coming from glutamine very low. the volumetric reaction rates of the ems were defined on the basis of the following assumptions: all reaction rates are specific, the metabolic reactions considered are all irreversible and therefore the respective reaction rates are nonnegative uncertainty in relation to kinetic constants and possible unknown saturation and inhibition effects. only the concentrations of glucose, glutamine and ammonia have a significant impact on the specific reaction kinetics. lactate never reaches inhibitory levels in our experiments. in the reaction rates of eq. the term in parenthesis represents the a priori knowledge concerning the kinetics of the particular reaction and ) whereas the ρi, μ and π terms represent the uncertainty concerning the reaction kinetics ) and are functions of three state variables ). with this particular formulation, the vector of known kinetic functions is given by: ϕ = t,     whereas the vector of unknown kinetics is given by: ρ = t = ρ     the last term in eq. is the control input vector u that in our case accounts for the volumetric feeding of glucose, fglc, and glutamine, fgln. identification of the em kinetics an important point for the identification of unknown flux functions is that the rank of k is thus the measurement of = > rank) is sufficient for the observability of the em kinetics. the other relevant point is the availability of sufficiently "rich" measurements to identify the "true" fluxes. preliminary simulation tests showed that, for the system structure of eq, the "true" fluxes can be identified under typical fed-batch conditions. experimental data of seven experiments were used for the identification of the em kinetics. data of experiments were used for parameter calibration whereas data of experiments were used for model validation. the concentrations in the state space vector were analyzed off-line according to methods published elsewhere. the neural network has three inputs, t and outputs as defined by eq.. the number of hidden nodes was tuned heuristically in the sense of minimizing the error of the validation data set. the best result was obtained with hidden nodes giving a total number of network parameters equal to the output scaling factors reflect the maximum kinetic rates and were defined as ρmax = t. the hybrid modeling results in terms of measured and predicted state variables are presented in figure for both training and validation data sets. examples of em kinetics identification are provided in figure the hybrid model was able to describe simultaneously all seven experiments with acceptable accuracy. in particular, the results with the validation data set strengthen the predictive potential of the model. metabolic interpretation in general, the metabolic activity of bhk cells during the cell growth phase is higher in the batch culture than in the fed-batch culture. the em fluxes in the fed-batch experiments seem to be much more controlled than in the batch experiments. the high levels of glucose and glutamine at the beginning of the batch culture are mostly directed toward the overflow metabolism, i.e., waste production of lactate and alanine. the fed-batch culture, which had lower glucose and glutamine concentrations during the cell growth phase, started with substantially lower consumption rates of these nutrients, while maintaining the flux to biomass. glucose is consumed for biomass synthesis and is metabolized through elementary modes em and em the consumption of this nutrient differs significantly from batch to fed-batch cultures. cells use glucose in a more efficient way in fed-batch than in batch cultures particularly at the beginning of the culture, since the glucose metabolized to lactate is much higher in the batch experiment. these results are in agreement with several published works for other mammalian cells, where it is reported that high glycolytic activity of animal cells results from high residual glucose. glutamine is consumed for product and biomass synthesis and is metabolized through elementary modes em and em the most energetic pathway involving glutamine, em is practically inactive during the cell growth phase of the batch culture, being glutamine preferentially converted into alanine. on the opposite, this elementary flux mode is an important pathway in the course of the fed-batch culture, representing 50% of the total glutamine consumption. the higher production rate of alanine at high glutamine levels is consistent with observations made by doverskog et al. and vriezen et al.. concerning the specific formation rate of glycoprotein, πigg, the product formation is consistently more stable in the fed-batch culture than in the batch culture. the product synthesis rate oscillates between and mg 10- cells h- in the former case whereas in the latter case the reaction is much lower in the beginning. this appears to be correlated with the overflow metabolism in the batch experiments, which seems to be detrimental for product synthesis. on-line culture optimisation the hybrid model was further used for on-line optimisation of a fed-batch bhk culture. the model parameters were re-tuned on-line using the data from the running cultivation. a batch training scheme was adopted, whereby the data of the running cultivation along with the data of historical experiments were used for model adjustment. some variables, namely glucose, glutamine, lactate and viable cells concentration, can be measured off-line, with the results available after a short period of time. therefore, at sampling times of to h, these measurements were stored in the training data set and then used for parameter identification using the same strategy previously described in the identification of the em kinetics section. the only difference was that the initial parameter values were those obtained in the off-line training procedure. the ammonia, alanine and igg concentrations were exclude form this model adjustment step because they were quantified only at the end of the experiment. after the parameters retuning step, the glucose and glutamine feeding rates were re-optimised in the sense of maximising the total amount of igg1-il produced at the end of the experiment. max⁡uj=cigg−il2vequation 0≤fglu≤fglu,max0≤fgln≤fgln,max−1∫t0tfriskdτ≤riskmax⁡      mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaafaqaaeqbcaaaaeaadawfqaqaaigbc2gatjabcggahjabciha4bwcbagaemydauhabeaaaoqaaiabdqeakjabg2da9iabdoeadnaabaaaleaacqwgjbqscqwgnbwzcqwghbwrcqghsislcqwgjbqscqwgmbatcqaiyagmaeqaaogaeiikagiaemidaq3aasbaasqaaiabdagambqabagccqggpaqkcqwgwbgvcqggoaakcqwg0baddawgaawcbagaemozaygabeaakiabcmcapaqaaaqaaiabbwgaljabbghaxjabbwha1jabbggahjabbsha0jabbmgapjabb+gavjabb6gaujabbccagmaabmaabagaegyoakdacagloagaayzkaaaabaaabagaegimaajaeyizimqaemoray0aasbaasqaaiabdeganjabdygasjabdwha1bqabagccqghkjyocqwggbgrdawgaawcbagaem4zacmaemibawmaemydaunaeiilawccbigae8xba0mae8xyaemae8heaghabeaaaoqaaaqaaiabicdawiabgsmijkabdaeagnaabaaaleaacqwgnbwzcqwfsbabcqwfubgbaeqaaogaeyizimqaemoray0aasbaasqaaiabdeganjab=xgasjab=5gaujabcycasiab=1gatjab=fgahjab=hha4bqabaaakeaaaeaacqggoaakcqwg0baddawgaawcbagaemozaygabeaakiabgkhitiabdsha0naabaaaleaacqaiwaamaeqaaogaeiykakyaawbaasqabeaacqghsislcqaixaqmaagcdawdxbqaaiabdkfasjabdmeajjabdofatjabduealnaadmaabaacbegae43yammaeiikagcccigae0hxdqnaeiykakcacaglbbgaayzxaagaemizaqmae0hxdqnaeyizimqaemouailaemysakkaem4uamlaem4sas0aasbaasqaaigbc2gatjabcggahjabciha4bqabaaabagaemidaq3aasbaawqaaiabicdawaqabaaaleaacqwg0baddawgaaadbagaemozaygabeaaa0gaey4kiipaaagccawljagaaczcamaabmaabagaegymaejaeg4mamdacagloagaayzkaaaaaa@aed9@ the optimisation is subjected to the constraints of the process dynamics, upper and lower bounds for the glucose and glutamine feeding rates, and the maximum risk of model unreliability, riskmax. the increase in volume is negligible in this problem, thus it was not considered in optimisation. the risk constraint states that the average risk must not exceed a given maximum level defined by the user. this restricts the feasibility domain to low risk regions and is essential for process optimisation supported by hybrid models since the black-box model predictions degrade in regions of the input space with sparse measurements. the micro-genetic algorithm coded by was used to solve optimisation. the population size and number of generations was and respectively. a maximum risk level of 35% was adopted during the on-line optimisation experiment. at each sampling time, the flow rates of the feeding pumps were updated according to the re-optimised feeding profiles of glucose and glutamine. medium optimisation higher productivities are likely if the initial medium composition is optimised. the optimisation results of initial glucose and glutamine concentrations along with the corresponding feeding strategies are shown in figure 8a. the medium should have initially low levels of glucose and glutamine. as for the feeding strategies, low levels of glucose and glutamine are optimal during the cell growth phase, whereas during the cell death phase, glutamine should slightly increase and glucose should decrease to a concentration of mm. a significant increase of product titre was predicted under such cultivation conditions. the kinetics of the elementary modes for the optimal strategy is presented in figure 8b. the optimal strategy can be interpreted by means of intracellular flux distribution. the results were not much different from the elementary flux distribution of the fed-batch culture. the main difference arises from the fact that the specific product formation increases during the cell death phase as a consequence of the increase of the glutamine concentration.  the current work presents a hybrid modeling approach for bioprocesses that integrates information concerning intracellular metabolic fluxes, bioreactor transport phenomena and measured data. first, elementary flux analysis is used to reduce the biosystem metabolic network into a set of macroscopic reactions relating extracellular components only. a bioreactor dynamical model is then established from material balance equations of the compounds which intervene in the final reaction set. nonparametric techniques are used for identification of the elementary modes kinetics from measured data. the method was successfully applied to a recombinant bhk-21a cell line producing a fusion glycoprotein. a significant result is the achievement of the flux distribution over the runtime of a bioprocess. the so obtained information can be used to identify conditions that favour product formation. a fed-batch bhk culture performed with on-line optimisation supported by the proposed methodology allowed a 10% increase in the final productivity. higher productivities are expected if starting nutrients concentrations are optimised. the developed tool promises to be advantageous for optimising the productivity of fed-batch biochemical processes since the transfer and adaptation to different cell lines is reasonably straightforward 
a dynamical hybrid system integrating material balance equations, ems reactions stoichiometry and kinetics was formulated. moreover, the method allows the identification of structural modifications in metabolism over batch time. in this work, a bioprocess control method is presented, which is designed on the basis of the metabolic network of the organism under consideration. the metabolic network of the system under study is decomposed into elementary modes, which are the simplest paths able to operate coherently in steady-state. the identified ems kinetics demonstrated typical glucose and glutamine metabolic responses during cell growth and igg1-il synthesis. a reduced reaction mechanism in the form of simplified reactions connecting substrates with end-products is obtained. the methodology was employed to analyse experimental data of recombinant baby hamster kidney cultures producing a recombinant fusion glycoprotein. ems kinetics were defined as the product of two terms: a mechanistic/empirical known term and an unknown term that must be identified from data, in a process optimisation perspective. finally, an online optimisation study was conducted in which the optimal feeding strategies of glucose and glutamine were calculated after re-estimation of model parameters at each sampling time. the progress in the "-omic" sciences has allowed a deeper knowledge on many biological systems with industrial interest. the main contribution of this work is a novel bioreactor optimal control method that uses detailed information concerning the metabolism of the underlying biological system. the bioprocess dynamics are formulated using hybrid rigorous/data driven systems and its inherent structure is defined by the metabolism elementary modes. this approach allows the quantification of fluxes carried by individual elementary modes which is of great help to identify dominant pathways as a function of environmental conditions.
 the progress in the "-omic" sciences has allowed a deeper knowledge on many biological systems with industrial interest. this knowledge is still rarely used for advanced bioprocess monitoring and control at the bioreactor level. in this work, a bioprocess control method is presented, which is designed on the basis of the metabolic network of the organism under consideration. the bioprocess dynamics are formulated using hybrid rigorous/data driven systems and its inherent structure is defined by the metabolism elementary modes.  the metabolic network of the system under study is decomposed into elementary modes, which are the simplest paths able to operate coherently in steady-state. a reduced reaction mechanism in the form of simplified reactions connecting substrates with end-products is obtained. a dynamical hybrid system integrating material balance equations, ems reactions stoichiometry and kinetics was formulated. ems kinetics were defined as the product of two terms: a mechanistic/empirical known term and an unknown term that must be identified from data, in a process optimisation perspective. this approach allows the quantification of fluxes carried by individual elementary modes which is of great help to identify dominant pathways as a function of environmental conditions. the methodology was employed to analyse experimental data of recombinant baby hamster kidney cultures producing a recombinant fusion glycoprotein. the identified ems kinetics demonstrated typical glucose and glutamine metabolic responses during cell growth and igg1-il synthesis. finally, an online optimisation study was conducted in which the optimal feeding strategies of glucose and glutamine were calculated after re-estimation of model parameters at each sampling time. an improvement in the final product concentration was obtained as a result of this online optimisation.  the main contribution of this work is a novel bioreactor optimal control method that uses detailed information concerning the metabolism of the underlying biological system. moreover, the method allows the identification of structural modifications in metabolism over batch time. 