Modelling and parameter estimation of gene expression and cell growth in batch cultures

Modelling and parameter estimation of gene expression and cell growth in batch cultures 1 R. Cubarsi Dept. Matem`atica Aplicada i Telem`atica Universitat Polit`ecnica de Catalunya, Barcelona, Spain J. L. Corchero, P. Vila and A. Villaverde Institut de Biologia Fonamental Universitat Aut`onoma de Barcelona, Bellaterra, Spain 1This work was presented at the 3ecm Third European Congress of Mathematics, held 10-14 July, 2000, Barcelona, Spain, and it has been supported by CICYT of Spain under grant No. BIO95-0801, and partially by Generalitat de Catalunya under Grant 1996XT-00030, and CUR (CIRIT) under grant 1995SGHR 00376. Modelling and parameter estimation of gene expression and cell growth in batch cultures R. Cubarsi1, J. L. Corchero2, P. Vila2and A. Villaverde2 1Dept. Matem`atica Aplicada i Telem`atica, Universitat Polit`ecn ica de Catalunya, Barcelona, Spain; 2Institut de Biologia Fonamental and Departament de Gen`etica i Microbiologia, Universitat Aut`onoma de Barcelona, Bellaterra, Spain Abstract Experimental procedure of CI857ts-controlled recombinant gene expression in bacterial batch cultures is mathematically modelled, and the corresponding minimum variance parameters are estimated from specific statistical or numerical methods, basically by using a global and recursive weighted least squares procedure under some constraints induced by the model. Moreover the numerical techniques proposed in this work act by accumulation of data coming from several runs of the experiment, so that more accuracy is obtained in the parameter estimation. In particular, for the production process, an extra-model parameter depending on an indicator vector is introduced for each run of the experiment in order to globalize the data. The analysis of obtained data leads to an integrated model for both cell growth and gene expression, which describes an asymmetric dynamics between culture growth and protein yield, and can serve to predict the maximal value of accumulated protein and the time required for it to be achieved at any stage of the preinducing cell growth. CORRESPONDING AUTHOR: Rafael Cubarsi Dept. Matem`atica Aplicada i Telem`atica Universitat Polit`ecnica de Catalunya, Campus Nord Jordi Girona, 1-3; E08034-Barcelona; Spain Phone: 34-3-401-5995, Fax: 34-3-401-5981 E-mail: [email protected]c.es 1991 MATHEMATICS SUBJECT CLASSIFICATION: 62, 65, 92. KEY WORDS: mathematical modelling, parameter estimation, constrained least squares 1 1 Introduction One of the mechanisms commonly used to estimulate the expression of the recombinant genes, and consequently, the production of the encoded proteins, is a rapid increase of temperature at which the cells are cultured. At the permissive temperature, 28◦C, there is no recombinant gene expression, but when the temperature is shifted to 42◦C, cells start the synthesis of the recombinant product while they are also growing in the culture. This temperature-mediated induction of gene expression, which is very convenient for industrial purposes, is achieved by the use of two kind of controllers of the gene expression, which are also introduced into the recombinant cells. They are a positive regulator, the lambda pLand/or pRpromoters, and a negative regulator, the repressor CI857, which is active below 32◦C but it becomes efficiently inactivated at 42◦C (Villaverde et al., 1993). We have developed mathematical procedures to analyze the performance of protein production in cultures of recombinant E. coli submitted to heat induction (Cubarsi et al. 1998). For this analysis, two types of mathematical procedures are required. The first type is composed of statistical and numerical methods for parameter and error estimation, fitting curves, etc. But when a function is approximated from a set of data, the problem of what kind of functions must be used always arises. In our case the function must be interpreted from a biological viewpoint, and it must also describe some biological properties of the experimental system. Thus, above techniques can be correctly used only if the biological system has been modeled, and the system properties to be quantified have been focussed. This is the other mathematical aspect of the work, in fact to be done previously to the first one. Under the assumption that cell growth is not significantly altered by the presence in the cell of the recombinant protein β-galactosidase, the dynamics of the gene expression is modeled in two steps, so that the culture growing model is combined with the gene expression model, a first order differential equation that describes the protein production in terms of cell growth, in order to explain the time evolution of recombinant protein yield along the induction phase. Then, the set of parameters for both models is estimated by using specific least squares techniques subject to constraints from the models, with statistical evaluation of error propagation. In order to minimize the errors, the numerical algorithms for parameter estimation take advantage of working with a batch culture procedure with multiple induction sequences of the culture, where data from different induction sequences of the same non-induced culture are processed all together, as a single experiment. However, data from several runs of the same experimental process can be pooled in a global data set only under some specific requirements. Thus, for the culture growing process this can be done if the time interval between two consecutive culture samples remains constant along all the process. For the gene expression process, some parameters are constant for all the sequences, namely the model parameters, and others are sequence dependent. In this case, for each run of the experiment, an extra-model parameter depending on an indicator vector is introduced, such as an initial condition for the production process. Hence the total number of estimation parameters is increased by the number of runs of the experiment. The resulting model for cell growth and synthesis of recombinant proteins reveals an asymmetric distribution of the biosynthetic potential of the cell, which is manifested by a preferential synthesis of recombinant proteins in aged, slowly growing cultures. In other words, both growing and production capabilities of culture cells are not equidistributed, since when the culture growth velocity decreases, the protein production velocity still increases up to its maximum value. Moreover, the proposed model also allows a prediction of the optimal optical density of a batch culture to be temperatureinduced in order to get a predetermined amount of recombinant protein with the minimum induction time. 2 2 Basic notation The mathematical notation used in the work is now introduced by describing the experimental procedure. An initial amount of culture y0, measured from its optical density at the wave length of 550 nm (OD), is growing at 28◦C (initial stage). At this stage the culture growth can be described in terms of a time parameter tby a function y(t), which is the solution of an autonomous first order differential equation generated by a phase velocity field vythat, as we shall see in the following section, will depend on two parameters A0and B0: dy(t) dt =vy(y(t), A0, B0) (1) Hence the solution of this equation may be explicitly expressed depending on the parameters, and the initial value y0=y(0), as y=y(t, y0, A0, B0) (2) At a time tfrom the beginning of the experiment, a sample of culture with OD y(t) is transferred to a prewarmed bath at 42◦C (induction stage). Then the growth rate of the culture changes and the production of β-galactosidase protein begins. The production is measured in enzymatic units per ml, referred as β-gal in this work. The growth process under the induction conditions has a similar behaviour as in the initial stage, but with other model parameter values, namely A1and B1. After a time xin the induction stage, the OD of culture y(t), that had been induced at a time tfrom the beginning of the experiment, varies along what we shall call the t-induction sequence according to a function yt(x), which is the solution of a differential equation, similar to Eq. 1, such as dyt(x) dx =vy(yt(x), A1, B1) (3) Hence, the solution can be written explicitly as a function of the model parameters, and the initial value y(t) = yt(0), in the form yt=yt(x, y(t), A1, B1) (4) On the other hand, for the gene expression process along the induction stage, that is the recombinant protein production, we assume that the protein has not any toxic effect neither for itself, nor for the culture, and depends only on the OD of the growing culture. Details of experimental procedure are given by Corchero et al. (1994), where the functional dependence of β-gal production in terms of the OD of the induced culture has been proved. Thus, along the t-induction sequence, if yt(x) is the OD of an induced culture, for a given induction time x, we can evaluate the production process by means of a function βt(x) = β(yt(x)), depending on whether the time evolution or the culture OD dependency of the product is emphasized. Thus βt(x) represents the yield of β-gal for the induction time xin the same induction sequence. The functional dependence β(yt) is studied in the following sections from an approximation model given by a first order differential equation depending also on two parameters c1and c2: dβ dyt =φ(yt, c1, c2) (5) where φis an arbitrary function of the specified arguments, whose solution can be written for each t-induction sequence from an initial condition ct 0, so that ct 0=β(0), in the form β=β(yt, ct 0, c1, c2) (6) 3 Notice that ct 0is a function of y(t), that can be implicitly given by β(y(t), ct 0, c1, c2) = 0, since at the begining of the induction sequence there is no protein yield. The sub-index referred to the t-induction sequence will be omitted when the context provides sufficient information. Finally, the production kinetics, the time evolution of the protein production, can be studied by composition of the differential processes expressed in Eq. 3 and Eq. 5, so that the corresponding generating field is dβt dx =dβ dyt dyt dx =φ(yt(x), c1, c2)vy(yt(x), A1, B1) (7) Then, the fuction that describes the protein production in terms of the induction time xcan be expressed by recursive substitution of Eq. 2 and Eq. 4 in Eq. 6. 3 Mathematical model In the working conditions, and for both experiments ( (a)E42, Table 4, with induction temperature of culture at 42◦C, and (b)E40, Table 5, with induction temperature at 40◦C) the growth rate of culture can be satisfactory described, before and during the induction phase, by using the equation of limited growth of population models (see e.g. Hirsch & Smale, 1974): dy(x) dx =y(x)(A+By(x)) (8) with Aand Barbitrary constants (A > 0 and B < 0). At low OD’s, the growing rate is nearly constant but, at the same time as the biomass is increasing, the exponential growth stops and the OD of the culture tends to the asymptotic value l=−A B(9) that depends on the growth conditions. On the other hand, according to Corchero et al. (1994), a growing culture which has been induced over a time xproduces an amount of β-gal β(x), that depends nearly in a quadratic way on the biomass y(x) of that culture. Thus, non constant rate of production, with reference to the culture growth, may be reflected along the induction stage, even though toxic effects of the recombinant protein are excluded from this work. Then the relationship between the amount of recombinant protein β-gal and the OD of culture can be written as follows, β(y(x)) = c0+c1y(x) + c2y(x)2(10) The corresponding differential behaviour, according to Eq. 5, will then have the form dβ dy =c1+ 2c2y(11) The meaning of this relationship, from a biological viewpoint, is now investigated by assuming that cell division is not influenced by β-galactosidase protein, and that c1and c2are parameters of the model. Along the induction stage the culture is growing according to Eq. 8, with a growth velocity vygiven by vy(y) = y(A+By) (12) Thus, the function vy(y) is a parabola with vertex at y=−A 2B, corresponding to an OD the half of the limit value given by Eq. 9, and also corresponding to the maximum growth velocity. Similarly, 4 for the production phase, a first approach could assume the same behaviour for the proteins as for the culture, that is, during the induction stage the increase of β-gal is proportional to the increase of biomass, if non-negative. Hence the production velocity of Eq. 7, namely vβ, could be expressed in this simple model as vβ(y) = k vy(y) (13) with ka positive constant. In fact, data from Tables 4 and 5 suggest that the increasing of β-gal is always associated with the increasing of biomass (Flickinger & Rouse, 1993). Furthermore, notice that the condition of increasing biomass leads to a working interval 0 ≤y≤lfor the culture. Nevertheless a more complex behaviour, consistent with Eq. 10, must be adopted, since the production protein rate with respect to the culture growth is not constant. More specifically, the OD of culture corresponding to the maximum production velocity, namely m, could be different from the OD of culture for the maximum growth velocity, y=l/2. Therefore, the general case of Eq. 7 must be considered, according to vβ(y) = dβ dy vy(y); 0 < y < l (14) taking into account that, if toxicity phenomena are not present in the induction stage the following inequality must be satisfied, in the working interval: dβ dy ≥0 (15) This situation is studied in the first order approximation given by Eq. 11, which enable us to explain in a simple way the asymmetry that the production velocity curve may have with respect to the culture growth velocity curve. 4 Asymmetry between production and growth In order to compare both velocity curves of Eq. 14, the function in the right hand side member of Eq. 11 will be denoted, according to Eq. 5, as φ(y) = c1+ 2c2y(16) Then Eq. 14 becomes vβ(y) = φ(y)vy(y) (17) The condition expressed by Eq. 15 implies φ(y)≥0 in the working interval, and then it is easy to deduce that: (a) In any case c1must be positive. (b) If c2= 0 both velocities are proportional and they have a common maximum at m=l/2. (c) If c2>0 the maximum production velocity is reached after the maximum growth velocity of the culture, and m > l/2. (d) If c2<0, since c1≥ −2c2yis fulfilled in the interval 0 ≤y≤l, and the maximum value of the right hand side member is held at y=l, then the following inequality must be satisfied: c1≥ −2c2l(18) Thus the maximum production velocity is reached before the maximum growth velocity of the culture, and m < l/2 is also held. The three functions involved in Eq. 17 are non-negative for values of the OD within the working interval, and in its bounds vβ(0) and vβ(l) are null. Furthermore, taking into account the polynomial form of vβ(y), it is easy to see that there is a single maximum on this interval. Thus, by assuming 5 Exp. l/2ǫ E42 1.38 ±0.05 0 E40 1.71 ±0.03 0.36 ±0.03 Table 1: Parameters describing the asymmetry between β-gal production and culture growth from Eq. 20. c26= 0, the relative position ǫof the abscissa mreferred to the value l/2, corresponding to the maximum growth velocity of culture vy(y), is introduced m=l/2 + ǫ(19) Then the abcissa of the maximum can be written, in terms of an auxiliary parameter α=l+c1 c2, as follows ǫ=sign(c2)1 |α|+q|α|2+ 3l2 l2 2(20) Above expression is also useful in order to see how far can the maximum moves around the central value y=l/2, being consistent with the condition of Eq. 15. Notice that |ǫ|is a decreasing monotonic function of |α|, and from Eq. 18 it is easy to see that |α| ≥ l. Therefore, from Eq. 20, by substitution of this minimum value of |α|, we get the admissible range |ǫ| ≤ l 6. Also, taking into account Eq. 19, we can conclude that our model enable us to explain a maximum production velocitiy in the following range of values 1 3l≤m≤2 3l(21) Notice that if ǫ > 0 the age for significant production is delayed towards high values of OD, while for low OD’s the production would be insignificant. If ǫ < 0 the behaviour is in the opposite way. 5 Production kinetics In this section, for a given t-induction sequence, we study the time evolution of the product content βt(x) = β(yt(x)) that is present in the culture at an age xof the induction stage. Remember that the culture with OD yt(x) has been induced at a time tfrom the beginning of the experiment. Following the proposed approximation, according to Eq. 8 and Eq. 11, we can write Eq. 7 as follows dβt dx = (c1+ 2c2yt(x))(A+Byt(x))yt(x) (22) where c1>0, A > 0 and B < 0. Then, when the induction time x→ ∞, the OD of culture tends to the limit lgiven by Eq. 9. Hence the function βt(x) tends to the asymptotic value lim x→∞ βt(x) = β(l) (23) That is, from a sufficient large interval of time, the product concentration becomes nearly stationary. Moreover, since the factor c1+ 2c2ytof Eq. 22 is non-negative in the working interval 0 ≤y≤l, this asymptotic value is the maximum yield that can be reached. Thus the function βt(x) does not have any maximum before reaching their asymptotic value, and, for any induction sequence, the kinetic of the product has a monotonic increasing curve along all the induction process. In order to obtain the function βt(x) we must take into account how the culture is growing before and during the induction stage, since the function yt(x) depends also on the OD of culture just at the beginning of the t-induction sequence, y(t) = yt(0). Thus we write, according to the solution of 6 Eq. 8, the relationship describing the biomass evolution of a culture that has been induced over a time x, yt(x) = A1y(t)eA1x A1+ (1 −eA1x)B1y(t)(24) The sub-index 1 is used to distinguish the induction stage, and the value y(t) represents the OD of the culture at the beginning of the induction sequence, according to y(t) = A0y0eA0t A0+ (1 −eA0t)B0y0 (25) In the latter equation y0is the initial amount of culture at the beginning of the experiment, and the sub-index 0 is used to distinguish the growing stage before the induction. Finally, the production of β-gal in terms of the induction time xis obtained from Eq. 10, also combined with Eq. 24 and Eq. 25, by assuming that in the beginning of the t-induction sequence there is not any significative amount of product, βt(x) = c1(yt(x)−y(t)) + c2(yt(x)2−y(t)2) (26) Sometimes the foregoing equation will be used in order to describe the time evolution of β-gal, and sometimes the production in terms of OD of the induced culture. Some consequences and features of these equations will be pointed out in the last section. 6 Culture growth parameters The algorithm to estimate the growth parameters Aand Bof Eq. 8 for the induction stage, as well as in the initial phase, is based on the time equidistribution of the OD samples, as it is shown in Table 4 and Table 5 for both experiments. By inverting that expression, a linear dependence between the inverse of the OD of two consecutive culture samples with an arbitrary time separation ∆x=x−x0is obtained, 1 y(x)=e−A(x−x0)1 y(x0)−B A(1 −e−A(x−x0)) (27) Hence, by defining a=e−A∆x b=−B A(1 −e−A∆x) zk=1 yk (28) and maintaining the time interval ∆xconstant for any couple of consecutive samples, Eq. 27 can be written in a simple and recursive way as follows: (ξk, ηk) = (zk−1, zk), ηk=aξk+b;k= 1,...,n−1 (29) Thus, if all the points (ξk, ηk) are graphically represented, a straight line is obtained. This equation will be used in order to compute the auxiliary parameters aand bfor the culture growth by means of a linear least squares approximation, and by taking into account the covariance matrix of errors. The procedure can be briefly described as follows. The OD measurements ykare obtained with independent errors ∆yk, with zero means and common variance σ2 OD. The errors of zkare evaluated according to the linear approximation from Eq. 28, ∆zk≃ −∆yk/y2 k, so that the accuracy of zkis given by the variance V(∆zk) = z4 kσ2 OD (30) 7 experiment E42 E40 stage initial induction initial induction a0.373 ±0.036 0.189 ±0.046 0.725 ±0.024 0.525 ±0.011 b0.279 ±0.029 0.295 ±0.024 0.089 ±0.046 0.139 ±0.005 cov(a, b) -0.0009 -0.0011 -0.0011 -0.00005 A0.987 ±0.095 1.665 ±0.244 0.643 ±0.066 1.290 ±0.042 B-0.438 ±0.062 -0.605 ±0.103 -0.208 ±0.1114 -0.377 ±0.017 cov(A, B) -0.0057 -0.0249 -0.0071 -0.0007 sOD 0.09 0.25 0.03 0.13 Table 2: Culture growth parameters, estimated errors and covariances of the estimates, from Eq. 38 and Eq. 8, with weighted RMS error sOD for culture OD from Eq. 42. The parameters aand bare related to the straight lines of Figures 1 and 2. Thus Eq. 29, a system of n−1 equations may be explicitly written with the corresponding errors as follows zk=a zk−1+b+δk;δk=a∆zk−1−∆zk;k= 1,...,n−1 (31) Hence the error associated with each equation depends on the unknown parameter aand, taking into account Eq. 30 and Eq. 31, the errors of Eq. 31 system are correlated according to the following covariance matrix of error vector ~ δ, E(~ δ~ δt) = σ2 ODVz=σ2 OD          z4 1+a2z4 0−az4 10... 0 −az4 1z4 2+a2z4 1−az4 2... . . . 0−az4 2... 0 . . ....−az4 n−2 0... 0−az4 n−2z4 n−1+a2z4 n−2          (32) where Edenotes the expectation and the symbol tmeans transpose. It is well known that the value σ2 OD is not necessary in order to estimate the parameters aand b. However, since the matrix Vzdepends on a, the estimation must be done iteratively, by revaluating Vzin each step, and by assuming the initial covariance matrix as the identity matrix. On the other hand, σ2 OD must be known in order to evaluate the covariance matrix of estimators V(a,b), then an unbiased estimator s2 OD of σ2 OD (Stuart & Ord, 1991, pp.723 and 737) is given by the following inner product, where ~η is the vector of experimental values, and ~η ∗is the predicted values vector, s2 OD =1 n−3(~η −a~η ∗)tV−1 z(~η −a~η ∗) (33) Note that this is equivalent to evaluate a weighted root mean square (RMS) error with respect to the inverted covariance matrix, over the number of degrees of freedom, n−3, of the problem. Finally, the covariance matrix of parameters Aand B,V(A,B), is obtained by error propagation approximation from the jacobian matrix J, and the covariance matrix V(a,b)(see e.g. Barlow, 1989) V(A,B)=J V(a,b)Jt;J=∂(A, B) ∂(a, b)(34) The estimates and corresponding errors of aand b, as well as of the parameters Aand B, are listed in Table 2. Also the weighted RMS error s2 OD is given for OD estimations. Note that least squares approximation provides us with unbiased estimators of the parameters and of sampling variances and covariances of the estimators, without assumptions concerning the forms of the error distribution (Stuart & Ord, 1991, p.716). This is only necessary when testing hypotheses about the 8