Dynamics of amino acid metabolism of primary human liver cells in 3D bioreactors

The kinetics of 18 amino acids, ammonia (NH3) and urea (UREA) in 18 liver cell bioreactor runs were analyzed and simulated by a two-compartment model consisting of a system of 42 differential equations. The model parameters, most of them representing enzymatic activities, were identified and their values discussed with respect to the different liver cell bioreactor performance levels. The nitrogen balance based model was used as a tool to quantify the variability of runs and to describe different kinetic patterns of the amino acid metabolism, in particular with respect to glutamate (GLU) and aspartate (ASP).


Introduction
Bioreactor technology for extracorporeal liver support using primary human liver cells has been developed within the last decade [1]. The design of a bioreactor for maintaining the hepatocyte's full functionality is of great importance. The used multi-compartment bioreactor consists of three interwoven, independent capillary membrane systems. Two of them (compartments 1 and 2, in the following aggregated to the 'perfusion compartment') provide decentralized plasma flow and the third one (compartment 3) provides oxygen supply to the cells, which are localized in the extracapillary space (compartment 4, 'liver cell compartment'). This threedimensional spatial structure represents an artificial equivalent of the hepatic vasculature at the lobular level. The bioreactor is integrated into a perfusion system that enables monitoring and control of system conditions (see Fig. 1). Previous studies have shown that primary liver cells reconstitute to tissue-like structures after inoculation into the bioreactor and that they maintain metabolic activities over several weeks [2][3][4]. This bioreactor therefore provides a valuable tool to analyze the dynamics and network structures of the physiological and molecular interactions of liver cells under standardized conditions that closely reflect the situation in the natural organ.
Recently, data from this bioreactor system has been analyzed statistically as well as by fuzzy cluster and rule based data mining and pattern recognition methods in order to identify early performance predictors for the bioreactor's long-term performance [5][6][7]. With respect to the prominent role of hepatocytes in the amino acid metabolism, the analysis and modeling of metabolic pathways of amino acids could provide important information on the functional state of hepatocytes cultured in vitro. In initial studies the role of amino acid regulation in extracorporeal liver support systems was analyzed [8]. More recently, three different modeling paradigms, i.e. correlation networks, Bayesian networks and systems of differential equations were applied to characterize the concentration profiles of six amino acids and related nitrogen-containing compounds in bioreactor cultures of primary human liver cells [6]. The present work reports on the model based analysis of the measured kinetics of 18 amino acids as well as ammonia (NH3) and urea (UREA) of 18 liver cell bioreactor runs.

Cells and bioreactors
Primary human liver cells were isolated according to a method described elsewhere [9] from human donor livers (n=18) that were not suitable for transplantation due to organ injury. After isolation, cells were cultured within bioreactors over 7-34 days. (This broad variability in duration is caused by their use for liver support from the seventh day onwards and by the performance. Low performance runs were finished earlier. High performance runs were continued for a longer period of time. However, only the first 6 days of the runs were analyzed in this study.) Concentrations of free amino acids in the culture perfusates were determined by an automated reversed phase high performance liquid chromatography system (RP-HPLC) with precolumn derivatization using the ophthaldialdehyde method [10]. NH3 and UREA concentrations were determined using routine clinical analyzers (Roche Diagnostics, Heidelberg, Germany).
With respect to the model based analysis, two compartments of the multi-compartment capillary membrane bioreactor system were considered: The 'liver cell compartment' with the volume V 2 =600 mL contains the liver cells in the inter-capillary space. The 'perfusion compartment' with the volume V 1 =900 mL supplies a stream through the inside of the capillaries. This perfusion stream carries the concentrations of the compounds that are supplied to or removed from the bioreactor. Decentralized mass exchange at low gradients is achieved by independent perfusion of two medium capillary systems, enabling different perfusion modes [10]. In this study, counter-directional medium flow was performed to facilitate rapid matter distribution. Due to the high flow rate of the perfusion stream (250 mL min À1 ), an almost ideal mixing within the perfusion compartment can be assumed (Fig. 1). To this perfusion compartment two time-variant inflow streams are added with the flow rates F A (t) and F B (t) as defined by Eq. (1). F A (t) follows a step function from F A1 =150 mL h À1 down to F A2 =50 mL h À1 switching at t A (=1 day). F B (t) switches from F B1 =0 up to F B2 =1 mL h À1 at the time t B (see Table 4). The outflow rate F 0 (t) to the waste equals the sum of both inflow rates. (This is realized by overflow keeping a constant pressure in the perfusion compartment).
The inflow rate F A (t) carries the 18 amino acids and NH3 at the concentrations c Ai (see Table 1). The inflow rate F B (t) carries only the amino acid aspartate (ASP) at the concentration c B15 =1,500 lmol L À1 , i.e. c Bi =0 for all i " 15. The time courses of the concentrations c 0i (t) in the outflow stream may be considered to describe the response of the medium to the inoculation of the bioreactor with cells. These concentrations are in steadystate at c 0i (0)=c Ai prior to inoculation. Due to the almost ideal mixing in the perfusion compartment, the concentrations c 0i (t) in the perfusion stream are the same as in the outflow stream.

Data
A data set with the elements c i,j,k (i=1,...,20; j=1,...,18; k=1,...,K) for 20 kinetic variables and 18 bioreactor runs at K time points t k was analyzed. The kinetics over the first 6 days of the bioreactor runs were analyzed here (t k =6 days). Concentrations of NH3 and UREA were measured daily (t k =0-6 days). Amino acid concentrations were determined up to the third day daily and every third day afterwards (t k =0, 1, 2, 3, 6 days). The samples for the measurement of the c i,j,k were taken from the waste of the liver cell bioreactor system (see Fig. 1) which was emptied daily after accumulation of the outflow stream over the period of 24 h. The model based analysis described focuses on the amino acid metabolism as quantified by the measured time series of the concentrations of 18 amino acids (i=1,..., 18; see Table 1) as well as of NH3 (i=19) and UREA (i=20).
Each run was labeled by L j 2 {'low', 'medium', 'high'} describing the performance with respect to the long-term maintenance of the functionality of the liver cells within the bioreactor. Seven runs were labeled 'high' (j=1,..., 7), seven runs 'medium' (j=8,..., 14) and four runs 'low' (j=15,..., 18). This performance labeling Fig. 1 Scheme of the liver cell bioreactor with the perfusion circuit, the two inflow streams and the outflow stream. Measured data were acquired from the waste was provided by an expert based on his assessment of altogether 99 variables that were measured to quantitatively characterize the system.

Differential equation system
The differential equation system (2.1-2.10) was developed to describe the measured kinetics of the 18 amino acids as well as of NH3 and UREA (Table 1). This model takes into account two compartments as described in Cells and bioreactors, the 'perfusion compartment' and the 'liver cell compartment'.
The Eq. (2.1) describe the dynamics of the components i in the perfusion compartment by four terms: The first and second term represent the fresh medium inflow with the volumetric rates F A (t) and F B (t), respectively, from the reservoir into the perfusion compartment with the volume V 1 . The reservoir concentrations c Ai of the components i are specified in Table 1. The volumetric rates F A (t) and F B (t) are specified by Eq. (1). The third term in Eq. (2.1) denotes the outflow from the perfusion compartment into the waste with the volumetric rate F 0 (t). The last term describes the diffusion between the perfusion and the liver cell compartment with its rate being proportional to the difference of the concentration c 0,i in the perfusion compartment and the concentration c i in the liver cell compartment. The parameter p 0 is proportional to the diffusion coefficient and assumed to be the same for all compounds to simplify parameter identification.
The Eqs. (2.2-2.9) formulate the changes of the concentrations c i in the liver cell compartment with the volume V 2 due to the exchange with the perfusion compartment (the first term on the right side of each differential equation) and the metabolic reactions (the following term(s) on the right side of each differential equation) as drawn schematically in Fig. 2. ð2:4Þ ð2:5Þ The model drawn in Fig. 2 and its Eqs. (2.1-2.10) were constructed based on the mean kinetics of 20 measured variables averaged over the seven high-performance runs H1-H7 as shown in Fig. 3. The kinetics of these seven individual runs are highly correlated with their correlation coefficients being greater than 0.8 for at least six of the seven high-performance runs and the 11 amino acids MET, SER, THR, ARG, HIS, GLY, TRP, LEU, VAL, PHE, and ILE [6]. The dynamics of these amino acids are characterized by monotonously, almost exponentially decreasing kinetics [6]. Their dynamics were therefore described by first order reactions as formulated by the second term in Eq. (2.2). The observed decreasing kinetics of amino acid concentrations is the net result of anabolic and catabolic reactions, where the catabolism of amino acids (which was considered here only) surpasses the amino acid synthesis and release by Table 1 Model variables c i and c 0 , i for the liver cell compartment and the perfusion compartment, respectively, their initial values c i(0) and c 0,i(0) before inoculation of the bioreactor with liver cells, the concentration c Ai in the inflow stream fed with the flow rate F A (t) according to Eq. (1) and the stoichiometric coefficients s i for nitrogen Name Model variables c i and c 0 , I can be interpreted mainly as amino acid uptake by transamination and oxidative deamination (see Table 2). In addition to the 11 amino acids with highly correlated and monotonously decreasing kinetics, the time courses of the concentrations of the amino acids ALA, LYS, and TYR were, for simplification, also modeled using Eq. (2.2). The kinetics of the amino acids ASP, ASN, GLU, and GLN were modeled in greater detail by Eqs. (2.3-2.6). Glutamate (GLU) is known to be the central compound of amino acid metabolism. Its kinetics were found in [6] to be highly correlated with ASP. ASP was fed by the inflow rate F B (t) (see Cells and bioreactors) to avoid or compensate the exhaustion of ASP which is essential for NH3 elimination via the urea cycle. The degradation of amino acids is interpreted in a simplified form by transamination to GLU for the amino acids LEU, HIS, ARG, VAL, TRP, PHE, ILE, ALA, TYR, and LYS (denoted by 'AAg' in Fig. 2). The amino acids MET, SER, GLY, and THR (denoted by 'AAn' in Fig. 2), however, are not transaminated forming GLU, but deaminated via different pathways forming NH3. The parameters with the indices 15-18, 21-23 as shown in Fig. 2 and Table 2 describe the metabolic fluxes between ASP, ASN, GLU, GLN, and NH3. They can be interpreted in a simplified way by reactions catalyzed, e.g., by the enzymes asparaginase, aspartate aminotransferase (AST), glutamate oxalacetate transaminase (GOT), glutamine synthetase (GS) and glutamate dehydrogenase (GLDH). In general, the biochemical reactions were assumed to be either of linear or bilinear nature (first or second order). In two cases, however, additional assumptions were required to obtain sufficient model fits: (1) For the aggregated modeling of the NH3 elimination via the urea cycle in the presence of ASP a non-linear Michaelis-Menten-type kinetics was used (see Eq. 2.8). This reflects the observation that the UREA formation rate does not significantly increase any more at strongly elevated ASP concentrations (see ASP and UREA in Fig. 3).
(2) Protein synthesis, i.e. anabolism, seems to be switched on or is strongly increased after the third day. This was modeled using Eq. (2.10). The total protein (PROT) was not measured, but a corresponding variable c 21 was hypothetically introduced by Eq. (2.9) to improve the model fit of the other 20 variables to the measured data for t>3 days.
The model parameters were identified by fitting the model to the measured data c i,j,k minimizing the scaled mean square error (mse) as defined by Eq.  Table 1 were averaged over 24 h (i.e. over the time interval of the accumulation of the bioreactor outflow in the waste where the samples for the measurements were taken from daily) and then compared with the measured data c i,j,k .
The differential equations were solved using a Runge-Kutta fourth order algorithm. The parameter identification by model fitting minimizing the mse was performed using a simplex search method. MATLAB tools (The MathWorks, Inc., Natick, MA, USA) were used for all calculations.

Results and discussions
The mean time profiles of the 20 measured variables each averaged over the high performance runs H1-H7 are shown in Fig. 3. The results of the model fitting to these averaged data are shown in Fig. 4 and the identified model parameters p 0 ,..., p 24 are listed in Table 2 (third column). The confidence intervals for the parameter values (fourth and fifth column in Table 2) were identified by repeated model fitting to the randomly disturbed time series data with means and standard deviations as shown in Fig. 3. The fit of the modeled kinetics to the measured data is acceptable for 18 of the 20 variables with the exception of lysine (LYS) and alanine (ALA). The initial increase of the LYS kinetics could be hypothetically explained by proteolytic activities. The final increase of ALA can be explained by the ALA aspartate transferase activity. The decreasing kinetics of 11 out of the 18 amino acids is explained by transamination and oxidative deamination forming GLU and NH3. The feeding of ASP after the third day results not only in an increase of the concentration of ASP but also of asparagine (ASN), GLU, and glutamine (GLN) due to the activities of the enzymes AST, GOT, GS, and asparaginase (alternatively or additionally the activity of asparagine synthetase could be included in the model). Table 3 lists the model parameter values p m with the respective model fitting error mse as identified by individual model fitting to the seven high performance runs. Figs. 5-7 show the results of the model fitting to the data for the individual high performance runs H1, H5, and H7, respectively. The kinetics of run H2 (not shown here) is very similar to the mean kinetics of the seven high performance runs shown in Fig. 4. The quasistationary concentrations of several amino acids, such as MET and LEU, are lower in run H5 (see Fig. 6) than in run H2 (compare Fig. 4). This results in a p 0 value, which is about 6.5· higher for run H5 than for run H2 (see Table 3). While the variability of the parameter values p 1 ,..., p 24 can be explained by differences in the liver cell material obtained from different individual donors (having different age, weight, liver damage, etc.), it is not entirely clear at this stage how to explain such different values of the parameter p 0 that represents diffusion and potentially other phenomena taking place across and beyond the membrane between the perfusion and the liver cell compartment. These issues are addressed in greater detail in [11] as well as in studies described in [12] to elucidate the dynamics of electrolyte distribution in the bioreactor without liver cells using tracer experiments and input/output systems analysis.
The model, which was developed using the kinetic patterns of high performance runs was also fitted to the data of the 11 runs with medium and low performance. The aim of these model studies is to elucidate causes of low liver cell bioreactor performance. The fit of the model with 25 parameters to 18 bioreactor runs results in an array of 25·18 parameters. Three selected parameters with their values as identified by fitting the model to the individual 18 runs are listed in Table 4. Considering a tolerable model fitting error mse of about 2% or less, the model fit proved satisfactory for six of the seven high performance runs (with the exception of run H7) and for two of the seven medium performance runs (M8 and M10). However, the model fit was not acceptable for the high performance run H7, the medium performance runs M9, and M11-M14 as well as all four low performance runs L15-L18. The main reason for the insufficient model fits is that the measured kinetics of leucine (LEU), isoleucine (ILE), valine (VAL), ALA, and other amino acids are not decreasing in these runs. This behavior is shown in Fig. 7 for run H7. In these cases some parameter values, e.g., for p 1 , were found to be zero (negative parameter values were set to zero; see Table 4). These zero parameter values cause that the medium and low performance runs displayed in Fig. 8 are artificially close to each other. Run H7 labeled as a high performance one appears to be a medium performance run according to the model simulation results (see Figs. 7 and 8).  3 Mean kinetics (± SD) of 20 measured variables (18 amino acids, NH3, UREA) averaged over the kinetics of the n=7 high performance runs H1-H7 (n=7) Figure 8 displays the identified model parameters of the 18 liver cell bioreactor runs after principal component analysis (PCA) of the parameter array for the first and second principal components in a biplot [13]. The first principal component represents 53% of the total parameter variance and the first and second principal components together represent almost 64% of this variance. Most of the parameters p 1 -p 14 (with the exception of p 2 , p 3 , p 9 ) representing transamination and oxidative deamination activities mainly constitute the first principal component of the parameter array (see Table 2: PC_1>0.2 for these parameters). These 11 parameters are correlated: The 55 correlation coefficients calculated between the 11 parameters  Fig. 3 (n=7). The simulated kinetics c 0i (t): (thin lines) are those obtained from the model (2). The simulated kinetics c 0i (t) were averaged according to Eq. (3) over the past 24 h (thick lines) in order to use them for model fitting to the measured data that were acquired from the waste (outflow accumulated over 24 h) Table 2 Parameters p m of the model (2) with interpretation, estimated values and confidence intervals [p mlow , p mhigh ]; enzyme activities: AST aspartate aminotransferase, GOT glutamate oxalacetate transaminase, GS glutamine synthetase, GLDH glutamate dehydrogenase; values and units of the parameters p m of the model (2) as identified by the model fit to the mean kinetics averaged over the seven high performance runs (Fig. 3) estimated over the 18 runs have values between 0.54 and 0.987 with the median of 0.87 and the mean of 0.84. The values of these 11 parameters are high for six of the seven high performance runs and low for run H7 as well as for the medium and low performance runs as shown in Table 4 for parameter p 1 . The value of the parameter p 1 is strongly correlated with the bioreactor performance (p-value 0.00017 by the two-sided t-test of p m values of high versus medium or low performance runs). The mean value of the parameter p 1 averaged over the high performance runs equals 2.84 (±1.77). The parameter is zero for all low performance runs. In Fig. 8 the medium and low performance runs are displayed in the left part describing low transamination and oxidative deamination activities whereas six of the seven high performance runs (with the exception of run H7) are displayed in the right part representing high transamination and oxidative deamination activities. Some parameters are highly correlated over the 18 runs with a correlation coefficient r of more than 0.98, e.g., the parameters p 11 and p 13 (representing the catabolism of methionine and glycine) as well as p 4 and p 7 . Also, the parameter p 1 is highly correlated (r>0.96) with the parameters p 4 and p 7 . The high correlation of the parameters p 1 , p 4 , and p 7 representing the transamination of LEU, ILE and VAL can be hypothetically explained by the fact that LEU, ILE and VAL are branched-chain amino acids which cannot be transaminated in hepatocytes. To transaminate these amino acids, the activity of non-parenchymal liver cells is required.
Among the high performance runs there are further individual differences that are displayed in Fig. 8 by a Fig. 7 Measured and simulated kinetics as in Fig. 5, but for run H7 Fig. 6 Measured and simulated kinetics as in Fig. 5, but for run H5 low second principal component PC_2 (e.g., high values of the parameters p 2 and p 22 for run H1) and a high PC_2 (e.g., low values of the parameters p 2 and p 22 for runs H2, H5, and H7). The parameters p 2 and p 22 represent the uptake of histidine and glutamine (see Table 2). As shown in Fig. 5 for run H1, the high parameter value p 2 results in a low stationary concentration of HIS and the high parameter value p 22 results in high NH3 and low GLU levels. These effects could be related to the amino acid transport system N (SN1) that mediates specifically the uptake of histidine and glutamine [14]. The parameter p 21 that represents the activity of GS was found to be very low for the high performance runs H2 and H5 (see Table 3). GS plays an important role in the spatial organization (zonation) of the liver and is exclusively expressed in pericentrally located hepatocytes [15]. The low GS activity in the high performance runs H2 and H5 could therefore be caused by a high proportion of periportal versus pericentral hepatocytes in the cell preparations for these two runs.

Conclusion
The kinetics of 18 amino acids and the related nitrogencontaining compounds NH3 and UREA in a primary human liver cell bioreactor were analyzed and modeled using a differential equation system. The model focuses on the kinetics of GLU and ASP as well as on the formation and elimination of NH3 and the synthesis of UREA. It describes the degradation of amino acids by transamination, oxidative deamination and other specific reactions. In addition, the activities of selected enzymes such as AST, GOT, GS, and GLDH as well as, in a more aggregated form, the activities of urea cycle enzymes are included. The differential equation system Table 4 Parameters for the 18 liver cell bioreactor runs: t B as used in Eq. (1); the star (*) for run 13 denotes that the feed F B (t) was switched on during the period from the second to the fourth day; the scaled mean square error (mse) and the parameters p 1 , p 2 , and p 22 as used in Eq. (2)   does not represent a fully mechanistic but rather a phenomenological model since essential metabolic activities had to be neglected because they cannot be identified based on the measured data. The differential equation system allows the analysis of a number of representative liver cell functions in terms of their kinetic behavior. The identification of the model parameters by fitting the model responses to the measured data was used to generate hypotheses about the causes of specific differences between the bioreactor runs. The model fits were found to be very satisfactory for eight high and medium performance runs. The model is however inadequate for low performance runs and with respect to the LYS kinetics also for high performance runs. This is probably caused by the neglection of proteases activities in the model that appear to be relevant for low and medium performance runs. Both, protein synthesis and degradation could not be modeled in detail due to the lack of representative protein measurements.
The applied model based analysis of data obtained from the bioreactor system can be used to quantitatively evaluate the functional state of liver cell cultures under high performance conditions intended for clinical application in extracorporeal liver support systems. The approach can also be used to study the effect of several exogenous factors, e.g., of hormones or drugs, on hepatocyte metabolism in vitro. The model based analysis methods applied here therefore provide suitable tools for in silico studies supplementing in vitro studies of hepatocyte functions in a systems biological way.