Improving the Estimation of Parameters of Penicillin Fermentation Models

Copyright © IFAC Computer Applications in Biotechnology, Osaka, Japan, 1998

IMPROVING THE ESTIMATION OF PARAMETERS OF PENICILLIN FERMENTATION MODELS M.T. Syddall' G.C. Paul' C.A. Kent

.,1

• Centre for Bioprocess Engineering, School of Chemical Engineering, The University of Birmingham, Edgbaston, Birmingham, B 15 2TT, United Kingdom

Abstract: Models for use in control and estimation applications should match the process as closely as possible. Fermentation process models are usually complex, containing many states and parameters. Obtaining accurate estimates of the parameters of such models is a costly and time-consuming process. Here we show a way of reducing the time and cost by designing optimal experiments for parameter identification. The method presented uses genetic algorithms to search for input profiles which optimise scalar functions of the Fisher information matrix, thus maximising the improvement in the parameter estimates that may be obtained from each experiment performed. The Penicillium chrysogenum penicillin-G fermentation, a secondary metabolite fermentation, is used as an example. Copyright © 1998 1FAC Keywords: Fermentation processes, Modelling, Parameter estimation, Genetic algorithms, Optimal experiment design

1. INTRODUCTION

same model was not practically identifiable, although more frequent sampling did help to improve confidence in the estimates obtained. The fact that algorithmic parameter estimation methods such as Marquardt methods do not necessarily lead to global optima was also mentioned.

The use of differential equation based physiological models in the design of optimal production strategies for penicillin fermentations (Lim et al. 1986, San and Stephanopoulos 1989) and in the design of advanced controllers for the fermentation (di Massimo et al. 1992, van Impe and Bastin 1995) has previously been reported in the literature.

These difficulties may be related to the geometry of the least squares error surface on which the parameter estimation is performed. Holmberg and Ranta (1982) showed that the niveau curves (contours of constant error) about the optimal parameter set were long and narrow.

The models on which these approaches were based contain parameters which must be estimated from experimental data, which can often cause difficulties. Nihtila and Virkunnen (1977) reported lack of confidence in parameter estimates obtained tuning bacterial models with data from batch fermentations. Holmberg (1982) demonstrated the theoretical identifiability of a model incorporating Michaelis-Menten kinetics, but went on to show that, given limited samples of noisy data, this

1

Series of experiments are often performed, generating incremental improvements in the quality of the parameter estimates. This is typically costly and time-consuming, particularly in the case of fermentation modelling, where a single experiment may take a week or more. Designing experiments so as to gain the maximum improvement in the parameter estimates could save time and

Author to whom correspondence should be addressed

23

money in the development of models for use in advanced fermentation control.

E({3)

= EI 13 =b + 1 + -({3 -

2

DEI

8{3' l3 =b ({3 - b) 2

b)' -8-E, I ({3 - b) 8{38(:J l3 =b

+ higher order

2. TUNING THE MODEL PARAMETERS

(5)

terms

Because E has a minimum at {3 = b ...

Throughout this work, the penicillin fermentation is considered as being described by a nonlinear differential equation based model of the following form.

8EI

=

0

8{3 l3=b and

x(t) = f(x(t), (3, u(t))

(1)

y(t) = g(x(t), (3, u(t))

(2)

is positive definite Neglecting terms above the second derivative, and substituting equation 4 into equation 5, we have

In the above, x(t) is a vector of model states, (3 is a set of time-invariant parameters, and u(t) is the vector of inputs to the model. The output of the model· is y(t); this second equation is used to relate measurements to the model states. (In our case, for simplicity, we assume y(t) = x(t).) The model structure used in this work is given in section 4.1.

2

1 E ~ -({3 - b)' -8-E, I ({3 - b) 2 8{38{3 l3 =b

which describes a hyperparaboloid. Surfaces of constant E are hence hypereIIipsoids. The second derivative of the error value, given in equation 3, with respect to the parameters, is given by the following equation (Eykhoff 1974).

In order to be able to make use of the model for practical purposes, it must first be tuned so as to most accurately represent the fermentation. This may be done using a least squares based optimisation routine, with an objective function of the following form. 1

(6)

(7)

n

2 2)m(ti) -

x(ti))'W(m(td - X(ti)) (3)

The second term vanishes close to the optimal parameter set, as liml3->b(m(td - x(td) = O.

In the above, E is the error value, m( ti) is a vector of measurement values at times ti, the summation is carried out for n sample times, and W is a weighting matrix. In our case, W is a diagonal matrix, with the maximum values of the measured states along the diagonal. The prime ' denotes vector or matrix transposition.

8x(t)/B{3 is given by the following equation, obtained by differentiating equation 1.

E=

i=!

d ()~g)

~

Bf Bx(t)

= Bx 8i3 +

8f(t) D{3

(8)

3. THE FISHER INFORMATION MATRIX AND OPTIMAL EXPERIMENT DESIGN

2.1 A geometrical interpretation of the errors The Fisher Information Matrix (FIM), forms the basis of several cri teria used in the design of optimal experiments for model identification (see Table 1.) For a derivation of the FIM, see Eykhoff (1974). In its discrete form, applicable for cases where measurements are taken at discrete sample intervals rather than continuously, the FIM may be defined as follows .

The error function E can be considered as a hypersurface given by

E = E({3) -

EI 13 =b

(4)

where b denotes the optimal parameter set, and hence EI 13 =b gives the minimum value for E. If we assume that the error surface is smooth and continuous with respect to the parameter values around the optimal parameter set, we can approximate the surface using a Taylor expansion around the optimum.

Comparing equations 7 and 9 shows that the the FIM is an approximation to the second derivative

24

Criterion A

simplified A C

o E

Formula min(tr(FIM 1)) max(tr(FIM)) min(tr(FIM)) max(dct(FIM)) max(Amin(FIM))

Interpretation minimise mean variance minimise mean variance minimises relative (mean) volume minimises ellipsoid volume minimises longest axis

min(cond(FIM) = ~::7:g~~i) spherical as possible Table 1. Criteria for optimal experiment design, derived from the Fisher Information Matrix (FIM) Amin and Amax are the minimum and maximum eigenvalues of the FIM. The above definitions are taken from Waiter and Pronzato (1990)

modified E

of the error surface, being the first term in the .

f.

expressIOn or

Here we are using GAs to search for fermentation inputs which optimise the experiment design according to the D and modified E criteria, two of the more commonly used criteria. (In our work, we have attempted to maximise the reciprocal of the condition number, rather than minimising the condition number itself-the two approaches are equivalent.) For use with the genetic algorithm, our input profile has been divided into a stepped input, with discrete portions having constant value . In this way, the input pattern is only determined by the values of the alleles being modified by the genetic algorithm. This stepped input profile is simple to specify using computer control, and may be applied by hand. The length of each 'step' in the input profile was chosen to be five hours as this was considered reasonable for manual implementation, should that be needed. The choice of a five hour interval also reduces the length of the strings and hence the size of the search space being searched by the genetic algorithm.

8 2 x(t,j fJ{3 u{3' .

The FIM has been used in the design of experimental conditions for estimating parameters of batch fermentations (Yoo et al. 1986) and fedbatch fermentations (Kalogerakis and Luus 1984) . Munack (1989) has shown that fed-batch fermentations are better, from the point of view of identifiability, than batch fermentations and his work was focussed on seeking out an input trajectory that made identification as robust as possible.

4. USING GENETIC ALGORITHMS TO SEARCH FOR OPTIMAL INPUTS The FIM depends on the model structure and parameter set used, on the input applied to the fermentation, and, in its discrete form , on the sampling interval used in obtaining measurements . Given a model structure and an estimate of its parameters, if the sampling rate is fixed, an optimal experiment may be designed by seeking out the input profile which maximises one of the design criteria given in Table 1.

4.1 Example - the Penicillin Fermentation

Munack (1989) stated that the gradient technique used in searching for an optimal input profile may have stopped in a suboptimal point, thereby finding a good local optimum, but not necessarily a global optimum . Genetic algorithms (GAs) have been shown to behave well on multi modal functions (Goldberg 1989), being less likely to become stuck in local optima than conventional optimisation techniques . GAs only need to calculate the objective function in the course of their search no use is made of derivatives, and G As may be used where the search surface is neither smooth nor differentiable.

The results of computer studies aimed at finding the best input for use in identifying the parameters of a penicillin fermentation model derived from that of Paul and Thomas (1996), given in Table 2, are presented here. This model is used to illustrate the application of the GA input optimisation technique to a complex fermentation model, for which conventional optimal control related methods would prove mathematically involved. The parameter set used in the model is given in Table 3.

In order to use genetic algorithms, the problem to be solved must first be encoded as a string, which the genetic algorithm acts on as it searches. In each generation, the algorithm evaluates the fitness of every string in the population, mates the strings according to their fitnesses (reproduction), exchanges information between pairs of strings randomly (crossover), and finally changes a small number of string elements with a low probability (mutation) .

In this work, the Genetic and Evolutionary Algorithm Toolbox for MATLAB (Pohlheim 1996) was used. An initial population size of 100, with 3 subpopulations was selected, running over a maximum of 100 generations, with the feasible input profiles bounded between 0 and 3.5/160 dm 3 /hT. The upper input limit was set to be the rate which would fill the working volume of the fermenter (7dm 3 in a lOdm 3 fermenter) over the course

25

Criterion

dXo J-LOXlS - - - -'nXo -dt Ko+S Kl+S dXl = J-LeXl S _ /1-0XlS + 'nXo _ J-L2 X 2P dt Ke + S Ko + S Kl + S dX 2 - - = /1-1 Vic - J-L2 X 2 + J-L3 X 2 dt dX3 r - - = /1-2 X 2P - J-La X 3 dt

(10)

o Modified E

(11)

OleJ-LeXoS Ke+S

- moXoS --Kl+S

(13)

(15)

dt

J-LLL(Xo + Xl) -dL = - -...:....;"---'--'----'-dt (KSL

(16)

+ L)(l + SI KsT)

dP

However, the graphs of the simulated fermentation results suggest that a practical fermentation carried out using the D-optimal input design could run into difficult.ies with low oxygen concentration from around 40 hours on . That the 'optimal' design gives rise to what may be a practically impossible situation could be because either the model does not describe the dissolved oxygen concentration or the maximum feed profile permitted is excessive .

(17)

dt dV -=F dt

(18)

Table 2. Model of Paul and Thomas (1996), simplified to increase simulation speed (equations 11 to 13), and extended to consider lactose present at the start of the fermentation (equations 15 and 16), X. - morphologically distinguished biomass fractions, S - glucose concentration, L - glucose concentration, P - penicillin concentration, V - broth volume, F - input feed rate, for clarity, dilution terms have been omitted Parameter J-Lu

J-Ll J-L2 J-L3 J-Le J-La J-Lp J-L! J-Lh Ko

KJ K2

Valu e 0.0333 0.0193 0.0535 2.026e-3 0.4092 0.0437 0.0287 0.241 0.0028 0.0352 0.0905 0.3017

Parameter Ke Kp Ki Olo Ole Olp

mo ml KSL KSI 1'1

E-optimal Input 3.00e19 1.30e19

criteria are compared with the values for a typical constant input feed profile. The two designed inputs produce better values for the criterion for which they were designed than does the constant feed profile (a greater D criterion value for the 0optimal design, and a greater E criterion value for the E-optimal design).

mlPcVicS K2+S

dL

D-optimal Input 5. 18e24 1.55e19

Table 4. Values of design criteria obtained using GA-designed input profiles

(12)

(14)

dt

Constant Input 1.64e4 1.30e20

The fact that both designed inputs produce 'better' values for the criterion for which they were not designed than for the constant input feed profile may be because the two experiment design criteria compared in this work are not entirely independent (det FIM = I1~ ~ ~::~: Geometrically, the D criterion is attempting to minimise the volume of the confidence ellipsoids, whilst the E criterion is attempting to improve the 'roundness' of the same ellipsoids .

,\ ).

Value 0.0788 0.0413 0.3656 1.85 1.83 0.85 0.0256 0.0242 0.0864 1.09ge-3 0.0122

6. CONCLUSIONS Genetic algorithms have been shown to be useful in designing optimal experiments for parameter estimation for complex, nonlinear fermentation models, for which optimal control based approaches to experiment design could prove involved.

Table 3. Table of parameters used in simulating the model of the fermentation (assumed to be 160 hours), starting from a typical initial volume (3.5dm 3 ).

In the future it is intended that the E-optimal experiment design will be implemented; as the current D-optimal design may encounter practical difficulties, this should not be implemented. There is also scope for work investigating the design of experiments using scaled parameter values, so as to obtain designs with equal percentage errors, as opposed to equal absolute error magnitudes, for all of the parameters.

5. RESULTS Both D and modified E optimal experiment designs improved over the course of the 100 generations used, giving quite distinct input profiles at the end. Figure 1 shows the D and modified E input designs produced, along with results of simulations performed using these input profiles. In Table 4, the designed inputs' values for both

Acknowledgement: The financial support of the 88SRC

is gratefully acknowledged.

26

D-optimal

.

.. . ... ..

0.02 LL

E-optimal

LL

0.01

X v.:l

Cl)

0 0 40 50 40 30 20 10 . . . . . . . . . . . 0 40 0 40

Cl..

80

120

160

. .........

80

120

160

.........

40

80

120

Cl)

Cl..

40

80

120

40

7 ··

. . . . . .

.

. ..

..

80

120

160

. ~

......

.......

40

80

120

160

40

80

120

160

20

5 .. · ···· .. · o~--~~----~----~----~

o~----~----~----~----~

o

0 0 50 40 30 20 10 0 0 40

0 0 10

160

5 ···· ..

". '

0.01

X v.:l

20 0 0 10

.......... .......

0.02

o

160

40

80

120

160

Time

Time

Fig. 1. Simulation results for D and modified E designs, F - input feed rate in dm 3 /hr, I:X. - total biomass conceIltration, S .- glucose concentration , P - penicillin concentration , all concentrations in glll/I Canadian Journal of Chemical Engineering 62 , 730- 737. Lim, H. C., Y . .I . Tayelo , .I . M. Modak and P. 130nte (1986). Computational algorithms for optimal feed rates for a class of fed-batch fermentation: Numerical results for penicillin and cell mass production. Biotechnology and Bioengineering 28, 1408-1420. Munack, Axel (1989). Optimal feeding strategy for identification of Monod-type models by fed-batch experiments. In: Computer Applications in Fermentation Technology. Elsevier. pp. 195- 203. Nihtila, Markku and .Iouko Virkunnen (1977). Practical identifiability of growth and substrate consumption models. Biotechnology and Bioengin eering 19, 1831- 1850. Paul, G. C. and C. R . Thomas (1996). A structured model for hyphal differentiation and penicillin production using Penicillium chrysogenum. Biotechnology and Bioengineering 51 , 558- 572. Pohlheim , H. (1996) . Genetic and evolutionary algorithm tool box for use with MATLAB - documentation . Technical report. Technical Uni-

7. REFERENCES di Massimo , Christine, Paul A. Lant, Aidan SaundeI's , Gary A. Montague, Ming T. Tham and A. .Iulian Morris (1992). Bioprocess applications of model-based estimation techniques. Journal of Ch emical Technology and Biotec/tTwlogy 53 , 265 -277. Eykhoff, Pieter (1974) . System Identification. Wiley. Goldberg, D. E. (1989) . Genetic Algorithms in Search, Optimisation, and Machine Learning. Addison-Wesley. New York. Holmberg, Andrea (1982) . On the practical identifiability of microbial growth models incorporating Michaelis-Menten type nonlinearities. Mathematical Biosciences 62, 23- 43 . Holrnberg, Andrea and .Iukka Ranta (1982) . Procedures for parameter and state estimation of microbial growth process models. A utomatica 18, 181 - 193. Kalogerakis , Nicolas and Rein Luus (1984). Sequential experimental design of dynamic systems through the use of information index .

27

versity Ilmenau. (http://www.systemtechnik.tuilmenau.de/-pohlheim/GA_Toolbox/index.html). San, Ka-Yiu and Gregory Stephanopoulos (1989). Optimization of fed-batch penicillin fermentation: A case of singular optimal control with state constraints. Biotechnology and Bioengineering 34, 72- 78. van Impe, J. F. and G. Bastin (1995). Optimal adaptive control of fed-batch fermentation processes. Control Engineering Practice 3,939- 954. Yoo, Young J., M. Marino-Galarraga, J. Hong and R.T. Hatch (1986). Experimental design for parameter estimation from batch culture. Biotechnology and Bioengineering 28 , 836841.

28