PIECEWISE-LINEAR APPROXIMATION OF NONLINEAR CELLULAR GROWTH

PIECEWISE-LINEAR APPROXIMATION OF NONLINEAR CELLULAR GROWTH Eduardo Mojica ∗,∗∗ Alain Gauthier ∗ Naly Rakoto-Ravalontsalama ∗∗ ∗

Universidad de los Andes, Cra 1 No.18A-10, Bogot´ a, COLOMBIA. e-mail: {ea.mojica70, agauthie}@uniandes.ed.co ∗∗ EMN,B.P.20722, 44307 NANTES, Cedex 03, FRANCE. e-mail: [email protected]

Abstract: In this paper, the piecewise-linear approximations of nonlinear cellular growth using orthonormal canonical piecewise linear functions are developed. The approximation method is presented and validated with satisfactory results. Future directions to hybrid control using this approximation are introduced. Copyright c 2007 IFAC  Keywords: Bioreactor, nonlinear dynamics, piecewise-linear (PWL) models

1. INTRODUCTION

The mammalian cells of Baby Hamster Kidney (BHK), which are used in the production of the vaccine against the foot-and-mouth disease when cultured on a well-defined medium comprising glucose, glutamine and other amino acids produce cells mass, monoclonal antibodies and waste metabolites such as lactate, alanine and ammonia. Mammalian cells display multiple steady states with widely varying concentrations of cell mass, desired products as well as waste metabolites (Namjoshi et al., 2003),(Mojica and Gauthier, 2005), (Mojica, 2005). This means that for identical input conditions to a fed-batch reactor, the outlet conditions change depending on how the culture is made fed-batch. Several nonlinear models have been developed for mammalian cells (Mojica, 2005). Here we address the peculiar features of mammalian cell growth in batch, fedbatch and continuous operating conditions and we construct a hybrid dynamical model capable to simulate the phenomenon.

In this paper, we develop a mathematical model using PWL (Piecewise-linear) approximation with Orthonormal High Level Canonical Piecewise Linear functions (Juli´ an et al., 2000). In particular, the class of all continuous PWL functions defined over a compact domain in Rm , partitioned with a simplicial boundary configuration is considered. The problem to find a PWL model given a nonlinear model has been previously treated in (Sontag, 1981), (Girard, 2002). In recent work, (Azuma et al., 2006) presents an approximation of nonlinear system using the Lebesgue piecewise affine approximation, it is applied to hybrid system modeling of gene regulatory networks. These methods are not satisfactory to our purpose of to implement a hybrid probing control (Mojica et al., 2007) based on value at vertices. This choice is motivated by several facts. First of all, this class of functions uniformly approximate any continuous nonlinear function defined over a compact domain Rn (see (Juli´an et al., 2000)). Moreover, the canonical expression introduced in (Juli´ an et al., 2000) uses the minimum and exact number of parameters, and it is the first PWL expres-

sion able to represent PWL mappings in arbitrary dimensional domains. As a consequence of this, an efficient characterization is obtained from the viewpoint of memory storage and numerical evaluation (Castro et al., 2001). Second, the approximation can be used in real implementations, the points taken from nonlinear model may be replaced for points taken from sensors or data directly from the process, namely, it addresses the problem of finding a piecewise linear (PWL) approximation of system where a reasonable number of measure samples of the vector field is available (regression set) (Storace and Feo, 2004). Third, in this alternative approach we deal with an approximation which is easier to handle than the nonlinear model, in fact, it can profit many tools developed for hybrid systems, e.g., the MLD model based approach (Bemporad and Morari, 1999). Finally, this CPWL is used in a model based control termed probing control in (Mojica et al., 2007), as a part of the future work to develop a hybrid probing control. The comparative analysis and error approximation between this new biological model and a nonlinear model developed first in (Mojica and Gauthier, 2005), (Mojica, 2005) is shown. This paper is organized as follows. In section II we present the process description and a brief description about a nonlinear model developed previously in (Mojica and Gauthier, 2005). In section IV we present the approximation structure proposed for biological models based on Canonical piecewise-linear (CPWL) approximation. In section V we give some simulation results for different dilution rates, and we verify the phenomenon of multiple states with the CPWL model through transient analysis and finally in section VI some conclusions are drawn.

2. PROCESS DESCRIPTION The mammalian cells BHK are used in the production of the vaccine against the foot-and-mouth disease, they are cultivated in bioreactor, which is located in the biotechnology laboratory LIMOR de Colombia S.A, in Bogot´ a, Colombia. It is a 2500 liters bioreactor, which operates in fed-batch mode, have temperature and agitation speed controlled. In this, some experiments were carried out to obtain information of several cellular cultures, which were used for parameters adjustment of the mathematical nonlinear model proposed in (Mojica and Gauthier, 2005), (Mojica, 2005).

2.1 Nonlinear Model The biological system class of nonlinear dynamics may be described using the following model

(Rewienski and White, 2003): · x(t) = f (x(t)) + B(x(t))u(t) y(t) = Cx(t)

(1)

where x(t) ∈ Rn is a vector of state at time t, f : Rn → Rn are a nonlinear vector-valued functions, B is a state-dependent n × m input matrix, u : R → Rm is an input signal, C is an n × k output matrix and y : R → Rk is the output signal. Our specific application, which was conducted with a developed model derived from the general state dynamical model described in (Bastin and Dochain, 1990) based on the mass balance, presents a description including two additional state variables; a substrate concentration, glutamine, which was included because like it was shown in (Namjoshi et al., 2003),(Mojica and Gauthier, 2005), (Mojica, 2005), glutamine apart from being the major source of energy also contributes to antibody production how we will show in simulation section, the transient analysis evidence this behavior. For further details about the nonlinear model refer to (Mojica, 2005). The nonlinear model used here is ⎧ ⎪ ⎪ ⎪ ⎪ · ⎪ ⎪ x1 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ · x2 ⎪ ⎪ ⎪ ⎪ · ⎪ ⎪ x3 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ · ⎪ ⎪ ⎪ ⎩ x4

 = μm

 x2 2 ki

x2 +x2 +ks

−

 = −k1 μm



= −k2 μm  = k3 μm

x2 2 ki x2 2 ki

x2 2 ki





x3 x3 +ka



x1 −

kdm kds (µm −kdm x4 )(kds +x2 )



x2 +x2 +ks



x2 +x2 +ks

x2 +x2 +ks







x3 x3 +ka x3 x3 +ka

x3 x3 +ka







x1 − x1 u.

x1 + (S1i − x2 )u x1 + (S2i − x3 )u

x1 − x4 u (2)

where x1 g/L is the cellular concentration, x2 g/L is the glucose concentration, μm 1/h is the maximum cellular growth rate, ks g/L is the glucose saturation parameter, x3 is the glutamine concentration, x4 is the lactate concentration, ka = 0.06 g/L is the glutamine saturation parameter, k1,2,3 are yield coefficients, 1.46, 0.9, and 0.06 respectively. S1i = 1.5, g/L is the initial glucose concentration in the feeding medium, S2i = 0.9 g/L is the initial glutamine concentration in the feeding medium, kdm = 0.012 is the maximum cellular death rate, and kds = 0.08 is the death saturation parameter.

3. BIOLOGICAL CPWL MODEL In this paper we deal with the PWL approximation technique developed in the past by (Sontag,

1981) and improved in the past few years by Juli´ an et al. (Juli´ an et al., 1999),(Juli´ an et al., 2000), and later it was extended by Storace and Defeo (Storace and Feo, 2004), to analyze the qualitative behaviors of dynamical systems characterized by PWL flows that approximate the flows of some academic examples. To guarantee that the dynamical behavior of the PWL-approximate flow will be faithful to the original system for any values of some significant parameters, we use the input variable u like a parameter, and simulate the PWL-flow. First, we give some definitions about CPWL (Canonical Piecewise-linear) based on the previous work in (Juli´ an et al., 1999),(Juli´ an et al., 2000).

characterized by the property that it produces a division of the domain into simplices. If S is defined as in (3), the space of all continuous PWL mappings defined over the domain S partitioned with a simplicial boundary configuration H is denoted by PWL[SH ]. At this point it is convenient to introduce the set VS of all simplex vertices contained in S. Depending on the number of coordinates different from zero, these vertices are organized in classes; for instance, all the vertices in S that have k(k ≤ m) coordinates different from zero are called class k vertices. The set VS plays an important role due to the fact that the function values at the vertices of VS constitute all the information that is needed to fully characterize any PWL function fp : S −→ R1 .

3.1 Canonical Piecewise Linear Functions

There are many possible choices for the PWL basis functions, each of which is made up of N (linearly independent) functions belonging to PWL N −dimensional linear space. We have that any basis can be expressed as a linear combination of the elements of the β basis, which is defined by recursively applying the following generating function γ(u, v) (Juli´ an et al., 1999),(Juli´ an et al., 2000):

The orthonormal definition of the PWL functions given before in (Juli´ an et al., 2000) is used in this paper to represent the nonlinear static mapping of cellular growth process in bioreator. With this in mind, a brief description of this representation and its most important characteristics are given next, for further details refer to (Juli´ an et al., 1999),(Juli´ an et al., 2000). For this approximation we have that given a nonlinear system (2) the following steps are necessary to obtain a CPWL approximation: Step 1. Given a compact hyperrectangle grid domain, order all vertices of that grid, Step 2. Group the vertices into simplicial cells, Step 3. Find CPWL basis functions, Step 4. Find the CPWL model approximation. In order to perform these steps we present the following definitions and methods: Definition 1. Let x0 , x1 , ..., xn be n + 1 points in the n−dimensional space. A simplex (or polytope) (x0 , x1 , ..., xn ) is defined by 

(x0 , x1 , ..., xn ) = x : x =

n i=0

μi xn

 where 0 ≤ μi ≤ 1, i ∈ {1, ..., n} and ni=0 μi =1.A simplex is said to be proper if and only if it can not be contained in a (n − 1) dimensional hyperplane. The representation proposed in (Juli´ an et al., 1999) requires the definition of a rectangular compact domain of the form S = {x ∈ Rm : 0 ≤ xi ≤ ni δk , i = 1, 2, ..., m} . (3) where δ is the grid size and ni ∈ Z+ , being Z+ the set of positive integers. Then this domain is subdivided using a simplicial boundary configuration H, a simplicial boundary configuration is

γ(x1 , x2 ) = {||x1 | + x2 | − ||x2 | − x1 | + + |x1 | + |x2 | − |x2 − x1 |}/4

(4)

Using (6), the functions γ 0 (x1 ) = x1 , γ 1 (x1 ) = γ(x1 , x1 ), γ 2 (x1 , x2 ) = γ(x1 , x2 ), and in general γ k (x1 , ..., xk ) = γ(x1 , γ k−1 (x2 , ..., xk ))

(5)

are defined. A distinctive property of (5) is that it possesses k nesting of absolute value functions, and accordingly it is said to have nesting level (n.l) equal to k. The first element of the basis is a constant term γ 0 (1). the remaining elements are formed by the composition of the function γ k (·, ..., ·), k = 1, 2, ..., m, with the linear functions πk,jk (x) =xk − jk δk k = 1, 2, ..., m, jk = 0, 1, ..., nk − 1. As a result, the basis can be expressed in vector form, ordered according to its n.l. as   T T T Λ = Λ0 , Λ1 , ..., Λm

(6)

where Λi is the vector containing the generating functions defined in [4] with i nesting levels. Accordingly, any fp ∈ PWL[SH ] can be written as

T

fp (x) = c Λ(x) (7)   T T T T where c = c0 , c1 , ..., cm , and every vector

Remark 3. It should n be noticed that a choice of δ such that δL  i=0 ei  = e, or equivalently, δ=

ci is a parameter vector associated with the vector function Λi . In order to obtain an orthonormal basis, it is necessary to define an inner product on PWL[SH ]. The new basis elements are linear combination of (6), that is Υ(x) =T Λ(x), and the matrix T may be obtained by two different methodologies. In this work we use the one built using the Granschmidt procedure as given in (Juli´ an et al., 2000). For finding the required approximation, we use a routine of [5] that finds a vector of parameters c that is the solution of the least square problem minx Ax − b2 , being A = ΥT (X), X the input matrix and b the output to be approximated in sparse format. In accordance with [5], the CPWL approximation of the nonlinear function g is defined as the function fp ∈PWL[SH ] satisfying fCPWL = Ac.

(8)

In this section, it is only necessary to know the values of f at the vertices because, as stated above, this is all the information needed to uniquely define a function belonging to PWL[SH ]. If the function f is continuous, the following results quantify the precision of the approximation. Lemma 1. If f is continuous in S, which is the union of nonoverlapping simplices, then (9)

where   max  =Δ ∈ S x0 , x1 ∈ Δ (fP W L − f (x)) .(10) max

In addition, if the function is assumed to be Lipschitz continuous in S, a useful relationship between the approximation error and the grid size δ can be obtained. Lemma 2. Let f : S → R1 be a function satisfying fP W L − f (x) ≤ L x1 − x0  ,

(11)

∀x1 , x0 ∈ S. Then, fP W L satisfies  n      fP W L − f (x) ≤ δL  ei  , ∀x ∈ S, (12)   i=0

(13)

guarantees that the approximation error is bounded as follows: fP W L − f (x) ≤ e.

3.3 Cellular Growth CPWL Model The cellular growth approximation scheme presented is based on the orthonormal canonical piecewise linear functions, therefore the idea is to approximate the right-hand side nonlinear function to a CPWL function for each differential equation. The PWL approximation based in a fourth order nonlinear model is obtained over the domain  S = y ∈ R3 : a1 ≤ x1 ≤ b1 , a2 ≤ x2 ≤ b2 , a3 ≤ x3 ≤ b3 , a4 ≤ x2 ≤ b4 , a5 ≤ u ≤ b5 .}

3.2 Analysis of CPWL Approximation: Error Estimation

fP W L − f (x) ≤ , ∀x ∈ S,

e n L  i=0 ei 

with a1 = 0.0 b1 = 4 a2 = 0.0 b2 = 4.5 a3 = 0.0 b3 = 2.6 a4 = 0.0 b4 = 1.8 a5 = 0.0 b5 = 0.3 The domain S is partitioned by performing m1 , m2 , m3 , m4 subdivisions along the state components x1 , x2 , x3 , x4 respectively, and m5 , subdivisions along the parameter component u, as the control variable. The coefficients were derived from a set of samples of f corresponding to a regular grid of n1 × n2 × n3 × n4 × n5 point over domain S. In particular, we will consider the following PWL approximation fCPWL of f : m1 = m2 = m3 = m4 = 6 , m5 = 3.n1 = n2 = 15, n3 = n4 = n5 = 10. Through several values of subdivisions we find that the least values are performing the systems well from both qualitative and quantitative points of view, we have to point out that the algorithm using orthonormal basis allow a lager number of inputs due to the ability of handling sparse matrices. Thus, it is possible to detect the simplices that contribute to the approximation, and it gives a measure of the nonlinearity of the function. It is possible to take smaller values of m in order to reduce the number of parameter but the quality of the approximation become worse.

x1−Cells[g/L] response for input case u=0.0

x2−Glucose[g/L] response for input case u=0.0

3

x1−Cells[g/L] response for input case u=0.031

5

2.5

x2−Glucose[g/L] response for input case u=0.031

2.5

5

2

4

3

1.5

3

2

1

1

0.5

0

0

Nonlinear model cpwl model

4

Nonlinear model cpwl model

2 1.5

2

1 Nonlinear model cpwl model

0.5 0

0

10

20 30 Time [h]

40

50

0

x3−Glutamine[g/L] response for input case u=0.0 3

10

20 30 Time [h]

40

50

x4−Lactate[g/L] response for input case u=0.0

0

10

20 30 Time [h]

40

50

0

0

x3−Glutamine[g/L] response for input case u=0.031

0.2

3

Nonlinear model cpwl model

2.5

1

Nonlinear model cpwl model

20 30 Time [h]

40

50

x4−Lactate[g/L] response for input case u=0.031 0.14

Nonlinear model cpwl model

2.5

0.12

0.15 2

10

0.1

2

0.08 1.5

0.1

1.5 0.06

1

1

0.04

0.05 Nonlinear model cpwl model

0.5 0

0

10

20 30 Time [h]

40

50

0

0

10

20 30 Time [h]

40

Fig. 1. Comparison of system response computed with 4th order nonlinear and cpwl models to the dilution u = 0.0 4. SIMULATION RESULTS

0.02

0

0

50

0

10

20 30 Time [h]

40

50

Nonlinear model cpwl model 0

10

20 30 Time [h]

40

Fig. 2. Comparison of system response computed with 4th order nonlinear and cpwl models to the dilution u = 0.031 /h Table 1. Start-up Conditions for Transient Analysis

4.1 Simulation Case: 4th Order Nonlinear and CPWL

S1 (Glucose g/L)

S2 (Glutamine g/L)

1

1.2

2.6

2

4.2

4.2

3

7.2

2.6

4

4.2

0.6

5

7.2

4.2

3 1 2 3 4 5

2.5 Cells Concentration

We focus on the simulation of a more realistic case, a fourth order nonlinear model. The simulation provides a plausible explanation of the different steady state and the state dependence of the matrix B. First, we simulate the CPWL model without input variable, dilution rate, in order to obtain a batch mode behavior. All the simulations of the 4th order model have as initial value x0 = (0.12, 4.2, 2.6, 0.0). Fig. 1 shows comparison between the nonlinear system and CPWL approximation obtained by numerically integration of the dynamical system having on the right-hand side of fCPWL for a dilution rate u = 0.

0.5

2 1.5 1 0.5 0

0

10

20

30 Time [h]

40

50

60

8 1 2 3 4 5

7 Glucose Concentration

In this case we can see the behavior of the four state variables. The CPWL model is qualitatively similar to the original system. For simulation purpose we increased the dilution rate to the optimal value, and the CPWL presents very good agreement. The dilution rate was then increased to reach a high dilution rate, where normally the system is not working, and the cellular concentration is too low. It is of particular interest to show that the model presented here is able to vindicate the multiple state stationaries. Fig. 2 shows the comparison between the nonlinear system and CPWL approximation for the case when the dilution rate is u = 0.031/h.

6 5 4 3 2 1 0

0

10

20

30 Time [h]

40

50

60

4.2 Transient Analysis of Cells Concentration

Fig. 3. Transients of Cells and Glucose Concentration for u=0.036 for different start-up conditions

Now, we consider the transient analysis for different initial conditions at the same dilution rate, we find several steady state. In order to further investigate this steady states, we established the

dilution rate, u = 0.036/h, and several initial conditions, with cellular and lactate initial conditions maintained unmodified, varying glucose and glutamine initial concentration, see Table 1.

50

We can observe clearly three steady states (see Fig. 3). Steady state 1 results from low levels of S1 as well as S2 maintained in a medium value, and in a low value, they are curves 1 and 4. Steady state 2 results from medium value of S1 and high value for S2 , as also S2 medium and S1 high, curves 2 and 3. Finally, the steady state 3 results from high levels of S1 and S2 , curve 5. It is remarkably the capability of the CPWL model to show the three steady states. We can see the direct relation between high concentration of glucose and glutamine in the initial conditions, produce high concentration of lactate which is appreciated in the cellular death (see Fig. 3).

5. CONCLUSION AND FUTURE WORK In this paper preliminary results on a Piecewiselinear approximation, based on orthonormal CPWL functions, of nonlinear cellular growth have been presented. It has also been proved that this structure allows us to approximate the dynamical behavior for different initial conditions in the transient analysis. Also, the nonlinear characteristic presented over cellular growth about three steady states is shown for the CPWL model capturing the partially substitutable and partially complementary nature of the two substrates, glucose and glutamine. This model will be useful to nonlinear control, in special to hybrid control. We can take a nonlinear system, after that the system is linearized, and finally the hybrid control is applied. One interesting point over this CPWL model is the fact that we used some points of the nonlinear system model, and in a real implementation, these points instead to take from nonlinear systems, they could be taken from the sensors or data from the process. Future work is focused on development an optimal control of this class of hybrid system to implement it on cellular cultures. This model although identified for the BHK experiments seems to be generally applicable to mammalian systems other than BHK, such as Hybridoma and Chinese Hamster Ovary (CHO) cells that display similar behavior (Namjoshi et al., 2003).

REFERENCES Azuma, S., J. Imura and T. Sugie (2006). Lebesgue piecewise affine approximation of nonlinear systems and its application to hybrid system modeling of biosystems. Proceedings of the 45th IEEE Conference on Decision and Control pp. 2128–2133. Bastin, G. and D. Dochain (1990). On-line estimation and adaptive control of bioreactor. Elsevier Science Publisher B.V. London.

Bemporad, A. and M. Morari (1999). Control of systems integrating logic, dynamics, and constraints. Automatica 35, 407–427. Castro, L., O. Agamennoni and C. D’Attellis (2001). Wiener-like modelling: a different approach. Proceedings of the 9th Mediterranean Conference on Control and Automation. Girard, A. (2002). Approximative solutions of odes using piecewise linear vector fields. Proceedings of the 5th International Workshop on Computer Algebra in Scientific Computing pp. 107–120. Juli´ an, P., A. Desages and B. D’amico (2000). Orthonormal high level canonical pwl functions with applications to model reduction. IEEE Trans. Circuits Syst. I 47, 702–712. Juli´ an, P., A. Desages and O. Agamennoni (1999). High-level canonical piecewise linear representation using a simplicial partition. IEEE Trans. Circuits Syst. I 46, 463–480. Mojica, E. (2005). Identificaci´ on y Control de un Proceso de Crecimiento Celular en Bioreactor. Master thesis. University de los Andes, http://biblioteca.uniandes.edu.co. Mojica, E., A. Gauthier and N. RakotoRavalontsalama (2007). Probing control for pwl approximation of nonlinear cellular growth. to be presented at the 13th IEEE Conference on Control Applications. Mojica, E. and A. Gauthier (2005). Desarrollo y simulaci´ on de un modelo para un cultivo celular en bioreactor. Proc. 2nd Congreso colombiano de Bioingenier´ıa e Ingenier´ıa Biom´edica. Namjoshi, A., W. Hu and D. Ramkrishna (2003). Unveiling steady-state multiplicity in hybridoma cultures: The cybernetic approach. Biotechnol. Bioeng 81, 80–91. Rewienski, M. and J. White (2003). A trajectory piecewise-linear approach to model order reduction and fast simulation of nonlinear circuits and micromachined devices. IEEE Trans. Computer-Aided Design 22, 155–170. Sontag, E. (1981). Nonlinear regulation: The piecewise linear approach. IEEE Trans. on Automatic Control 26, 346–358. Storace, M. and O. De Feo (2004). Piecewiselinear approximation of nonlinear dynamical systems. IEEE Trans. Circuits Syst. I 22, 830–842.