METHODS OF DESIGN OF OPTIMAL EXPERIMENTS WITH
APPLICATION TO PARAMETER ESTIMATION IN
ENZYME CATALYTIC PROCESSES

We consider the problem of identifying enzyme operating stability which is very important e.g. in industrial production of enzymes or in understanding enzymatic pathways in a living system. Traditionally enzyme stability is determined empirically by time-consuming experiments: under constant reaction conditions (temperatures) the half-life of the catalyst is measured [1]. This procedure is repeated under different temperatures. Another approach for faster identification of enzyme operating stability is suggested in [2]. In an instationary operated reactor the reaction temperature is increased linearly and parameters of a kinetic model based on Arrhenius equations are estimated. Then the stability properties of enzymes are estimated as functions of the kinetic parameters. Experience shows however that such experiments do not deliver enough data to identify the kinetic parameters reliably.

Optimal experiments (temperature profiles) are necessary that allow the estimation of the kinetic parameters of enzyme catalytic processes and thus stability properties. This paper focuses on numerical methods for parameter estimation and design of optimal experiments.

The paper is organized as follows. Firstly, we describe parameter estimation problems in ordinary differential equations as well as specifying the model of enzyme reaction under consideration. Secondly, the method used for parameter estimation is described. It follows a very successful solution approach, namely the boundary value problem (BVP) approach, which combines infeasible point optimization algorithms like the generalized Gauss-Newton method with BVP techniques like multiple shooting to discretize the dynamic system. The discretized equations are then treated as nonlinear equality constraints. We solve the resulting finite-dimensional, possibly large-scale, nonlinear constrained least squares problem with the generalized Gauss-Newton method. Under certain regularity assumptions the method shows the final linear rate of local convergence. To ensure global convergence, a newly developed efficient step-size strategy is used. Special attention is paid to the discussion of why the generalized Gauss-Newton method is particularly appropriate for parameter estimation and does not converge to solutions with large residuals. The quality of parameter estimates is defined by statistical analysis, and the results of parameter estimation for an immobilized enzyme system finish this part of the paper. Thirdly, theoretical justification methods for the design of optimal experiments as well as numerical and experimental results are presented. Optimal experiments are designed for the nominal values of parameters which are known only to lie in a confidence region. Fourthly, the design of experiments that are less sensitive to parameter uncertainty is described. The paper is finished by conclusions.

which are subject to measurement errors ε_ij. Here ptrue are the ''true'' values of the parameters. Note, that several model quantities can be measured at a time t_j. According to the common approach, in order to determine the unknown parameters an optimization problem is solved. The possible constraints of this problem describe the specifics of the model (constraints on the initial and terminal states, constraints on parameters, etc.) and can be formally written as constraints at times

A typical solution approach to parameter estimation which is found very often in practice is the initial value or single shooting approach: the ODE system is repeatedly solved as an initial value problem, and unknown parameters including possibly initial values are iteratively improved by some optimization procedure.

In contrast to that, our numerical solution of the parameter estimation problem is based on the Boundary Value Problem (BVP) approach going back to [4]. The basic idea consists in parameterizing the dynamic equations (initial or boundary value problem) like a boundary value problem (e.g., by multiple shooting) and then performing simultaneously (in one iteration loop) the minimization of the cost function and the fulfilment of the constraints given by the discretized boundary value problem. It has been shown [4,5], that BVP methods (based on multiple shooting or collocation) are much more stable and efficient than the single shooting approach when solving parameter estimation problems.

Multiple shooting

The scheme of multiple shooting consists of the following. First one chooses a suitable grid of multiple shooting nodes τ_j

 

covering the interval where measurements are given.

At each grid point the values of the state variables s_j are chosen as additional unknowns and m ODE initial value problems

 

are solved on each subinterval I_j:=[τ_j, τ_j+1] to yield a solution y(t; s_j, p, q, u) for . Solutions of dynamic systems, generated by this procedure, are usually not continuous at τ_j. This has to be enforced by additional matching conditions. Inserting the computed values y(t, s_j, p, q, u), into problem (1) one obtains a constrained optimization problem in the variables (s, p) := (s₀,...,s_m, p)

 

Multiple shooting possesses several advantages:

1. It is possible to include a priori information about the state variables, e.g. from the measurements or from expert knowledge, by a proper choice of initial guesses for the additional variables s_j. Thus, it can be ensured that the initial solutions y(t;s_j, p, u) remain close to the observed data. It can be shown that this damps the influence of poor parameter guesses.

2. The adequate choice of initial guesses for the state variables (and the application of a Gauss-Newton method for the solution of the constrained least squares problem) typically avoids convergence to local minima with large residuals.

3. The scheme is numerically stable. The splitting of the integration interval limits error propagation and makes it possible to solve parameter estimation problems even for unstable or chaotic systems.

4. The BVP discretization induces a very specific structure in the problem equations, which can be exploited in particular for parallelization.

Gauss-Newton method

Parametrization of the dynamics yields a finite dimensional, possibly large-scale, nonlinear constrained least squares problem which can be formally written as

 

 

Note, that the equalities F₂(χ) = 0 include the matching conditions induced by multiple shooting. We assume that the functions F_i : , i = 1, 2, are twice continuously differentiable. The number of variables in problem (5) is equal to number of differential equations multiplied by the number of multiple shooting nodes plus the number of parameters. To solve problem (5) we use a generalized Gauss-Newton method according to which a new iterate is (basically) generated by

 

where the increment Δχ is the solution of the following linear constraint l₂ problem at χ = χ^k

 Here, J₁(χ) and J₂(χ) denote the Jacobians of F₁(χ) and F₂(χ) respectively.

then a linearized problem (7) has a unique solution Δχ^k and a unique Lagrange vector λ^k satisfying the following optimality conditions

 

Using (8) one can easily show that Δχ^k can be formally written with the help of a solution operator J⁺

 

The solution operator J⁺ is a generalized inverse, that is it satisfies J⁺J J⁺ = I, and is explicitly given by

 

Choosing the step lenght t^k by means of classical line search methods based on the exact penalty function

 

with sufficiently large weights α_i < 0, i = 1,..., m₂, ensures global convergence. However, it is well known that already in mildly ill-conditioned problems such a step-size strategy may be very inefficient since it may produce very small step-sizes. Therefore we use the "restrictive monotonicity test", see [5,6], that has proved to be very effective in practical applications.

Note, that the generalized Gauss-Newton method (6) - (7) has several advantages. First it does not use second order derivative information and the local linearized problems are linear constrained least squares problems. Under certain regularity assumptions at the solution, the method shows a good linear rate of local convergence. There are problems, however, for which the Gauss-Newton method may have a rather slow local convergence rate or may even fail. The reason is that the linearized model (7), which forms the basis of the Gauss-Newton method, is an inadequate representation of such nonlinear problems, since the second-order information cannot be ignored. These problems are called problems with large residuals.

Using SQP-type methods for the nonlinear constrained l₂ parameter estimation problem one could force convergence to a solution even in such a case. However, such solutions are undesirable in a certain sense. We can show that a solution of this type, even if it is a strict minimum of the nonlinear constrained l₂ problem (5), cannot be expected to be a continuous deformation of the "true" parameter values under perturbations caused by the measurement errors. Thus, slow local convergence of the full-step Gauss-Newton indicates deficiencies in the model or lack of data and can be considered as an advantage of the method. For a detailed analysis see [5].

 

Solving the linear l₂ problem

At each iteration of a Gauss-Newton method a linear least squares problem (7) has to be solved.

First, we make use of the special block structure of the Jacobian J(χ) which is induced by multiple shooting:

 

Every block column corresponds to the derivatives with respect to the discretization variables and parameters in one subinterval. The block rows with G-matrices are the derivatives of the continuity conditions

 

The block rows with D-matrices correspond to the derivatives of the functions F_i of the cost functional and the constraints of the nonlinear problem (5) excluding the continuity conditions. We use a fast, stable and efficient structure exploiting decomposition (see [4, 5]) to reduce the large linear least squares problem to a linear least squares problem with smaller dimension.

The number of variables in the resulting problem is equal to the number of parameters plus the number of differential equations. This so-called condensed problem may be solved using the methods described in [7].

 

Statistical sensitivity analysis for the estimates

An important question in parameter estimation is how good the computed estimates are. The answer is provided by sensitivity analysis. If the experimental data is normally distributed then the estimated solution χ^* of the parameter estimation problem is also a random variable which is normally distributed in the first order χ^* ~ N(χ^true, C) with the (unknown) true value χ^true as expected value and the variance-covariance matrix C given by

 

The variance-covariance matrix describes the confidence ellipsoid which is an approximation of the nonlinear confidence region of the estimated variables.
The 100(1- a)% linearized confidence ellipsoid can be described by (see [5])

 

Here, the probability factor where is the quantile of the χ² distribution. The linearized confidence ellipsoid G_L(α,χ^true,q,u) is contained exactly in a box determined by the confidence intervals

 

where . Here, C_ii denotes the diagonal elements of the covariance matrix C. The values are known as standard deviations of the variables χ_i.

Results of parameter estimation for Candida antarctica lipase on ionic resin ("Novozym")

Applying the methods of parameter estimation described earlier, "Novozym" shows that the data derived from the initial experiment does not deliver enough information to identify all parameters, see Table 1. The settings of the initial experiment are defined by the temperature shown in Fig. 2. The parameter estimation results show that the information received in the experiment is not at all enough to estimate parameters reliably. Thus we need to design additional experiments in order to provide sufficiently good data. The theoretical justification and methods for design of optimal experiments are briefly described in the next section.

 

 

---

Optimum Experimental Design: Problem Statement and Numerical Methods

In this section we consider the problem of nonlinear optimum experimental design and discuss a solution approach to this problem.

 

Experimental design optimization problem

Since the experimental data is randomly distributed, the estimated parameters are also random variables. Data evaluated under different experimental conditions leads to estimations of parameters with good or poor confidence regions depending on the experimental conditions. We would like to find those experiments that result in the best statistical quality for the estimated parameters and at the same time satisfy additional constraints e.g. experimental costs, safety, feasibility of experiments, validity of the model etc.

The aim of optimum experimental design is to construct N_ex experiments by choosing appropriate experimental variables
q = (q₁,...,q_Nex), and experimental controls u = (u₁,..., u_Nex) in order to maximize the statistical reliability of the unknown variables under estimation. For this purpose we optimize a design criterion depending on the variance-covariance matrix C (9).

The optimum experimental design approach leads to an optimal control problem in ODE systems of the following form

 

The constraint (11) describes control and path constraints for each experiment of N_ex experiments. Free variables of the optimization problem are the control profiles u_i(t) (e.g. temperature profiles of cooling/heating) and the time-independent control variables q_i (e.g. initial concentrations, properties of the experimental device) for all experiments.

Since the "size" of a confidence region is described by variance-covariance matrix C any suitable function of matrix C has to be taken as a cost functional. Typical choices are

 

 

Numerical methods

The experimental design optimization problem is a nonlinear constrained optimal control problem. The main difficulty lies in the non-standard objective function which is nonseparable and implicitly defined on the sensitivities of the underlying parameter estimation problem, i.e. on the derivatives of the solution of the ODE system with respect to the parameters and initial values.

The numerical methods are based on the direct approach, according to which the control functions are parametrized on an appropriate grid by support functions, the solution of the ODE systems and the state constraints are discretized. As a result we obtain a finite-dimensional constrained nonlinear optimization problem which is solved by an SQP method. The main effort for the solution of the optimization problem by the SQP method is spent on the calculation of the values of the objective function and the constraints as well as its gradients.

Efficient methods for derivative computations combining internal numerical differentiation [5] of the integration method and automatic differentiation of the model functions [8] have been developed, see [9,10,11]. For more detailed discussion of the numerical methods for nonlinear optimum experimental design see [11,12].

 

Results of parameter estimation for Candida antarctica lipase on ionic resin ("Novozym") with data from additional experiments

Using the methods of optimum experimental design we have optimized 5 additional temperature profiles. The following set-up has been chosen for the experiments: the duration of each experiment is 20 h; measurements are taken every half-hour; the temperature profile is parametrized as follows

 

 

Additionally, there are the following restrictions on the design parameters T_i and slope_i.

 

 

Results of parameter estimation for "Novozym'' with data from the initial and additional experiments are presented in Tables 2 and 3. Table 3 shows the estimated values of half-life over an interesting temperature range. The optimal temperature profiles are presented in Fig. 3. Figure 4 presents the fits before and after parameter estimation. The results of parameter estimation show that now we can estimate all half-lives with standard deviation below 30 %.

 

 

 

Figure 5. Sequential approach to experiment design.