Multi-Objective Optimization of Experiments Using Curvature and Fisher Information Matrix

Abstract

The bottleneck in creating dynamic models of biological networks and processes often lies in estimating unknown kinetic model parameters from experimental data. In this regard, experimental conditions have a strong influence on parameter identifiability and should therefore be optimized to give the maximum information for parameter estimation. Existing model-based design of experiment (MBDOE) methods commonly rely on the Fisher information matrix (FIM) for defining a metric of data informativeness. When the model behavior is highly nonlinear, FIM-based criteria may lead to suboptimal designs, as the FIM only accounts for the linear variation in the model outputs with respect to the parameters. In this work, we developed a multi-objective optimization (MOO) MBDOE, for which the model nonlinearity was taken into consideration through the use of curvature. The proposed MOO MBDOE involved maximizing data informativeness using a FIM-based metric and at the same time minimizing the model curvature. We demonstrated the advantages of the MOO MBDOE over existing FIM-based and other curvature-based MBDOEs in an application to the kinetic modeling of fed-batch fermentation of baker’s yeast.

1. Introduction

Dynamic models of biological networks and processes are often created to gain a better understanding of the system behavior. The creation of dynamic biological models requires the values of kinetic parameters, many of which are system-specific and typically not known a priori. These parameters are commonly estimated by calibrating model simulations to the available experimental data. Such parameter fitting is known to be challenging, as there often exist multiple parameter combinations that fit the available data equally well; that is, the model parameters are not identifiable [1,2,3,4,5]. While there exist a number of reasons for such lack of parameter identifiability, experimental conditions have a strong influence on this issue and thus should be carefully designed. In addition, biological experiments and data collection are often costly and time-consuming, further motivating the need for well-planned experiments that would give the maximum information given finite resources.

Model-based design of experiments (MBDOEs) offer a means for integrating dynamic modeling with experimental efforts, as illustrated by the iterative procedure in Figure 1. The role of the model here is to capture the knowledge and information about the system up to a given iteration. By using MBDOEs, one could harness this knowledge to guide experiments in the next iteration. MBDOE techniques have been used extensively for chemical process modeling [6], and more recently, they have been applied to the modeling of cellular processes [7,8]. For the purpose of parameter estimation, experiments are generally designed to improve the precision of the estimated parameters. In this regard, the Fisher information matrix (FIM), whose inverse provides an estimate of the lower bound of parameter variance-covariance by the Cramér-Rao inequality [9], has been commonly used to define the objective function in the optimal experimental design (see [8] and references therein). Since the turn of the century, FIM-based MBDOE methods have had newfound applications in the emerging area of systems biology [10,11,12,13,14,15,16]. Besides FIM, Bayesian approaches have also been used for MBDOEs, where given the prior distribution of the model parameters, experiments are designed to minimize the posterior parameter variance [8]. Bayesian MBDOE strategies have been applied to the modeling of biological networks for reducing parametric uncertainty [17,18,19]. While our work is concerned with MBDOEs for the purpose of parameter estimation, MBDOE strategies have also been developed and applied for discriminating between biological model structures [20,21,22,23,24] and reducing cellular process output uncertainty [19,25,26].

In this work, we focused on FIM-based MBDOEs for parameter estimation. The FIM relies on a linear approximation of the model behavior as a function of the parameters. More precisely, the FIM is computed as a function of the first-order parametric sensitivity coefficients (Jacobian matrix) of model outputs. For systems with a high degree of nonlinearity, the optimal experimental design using the FIM may perform poorly [27]. For this reason, Bates and Watts proposed a MBDOE based on minimizing model curvature by using the second-order parametric sensitivities (Hessian matrix) [28]. Hamilton and Watts further introduced a design criterion, called Q-optimality, based on a quadratic approximation of the volume of the parameter confidence region [29]. More recently, Benabbas et al., proposed two curvature-based MBDOEs [30]. In one design, the authors used a minimization of the root mean square (RMS) of the Hessian matrix, while in another design, they employed a constrained optimization guaranteeing the RMS to be lower than a given level. While the second strategy using a curvature threshold was demonstrated to give more informative experiments, how to set the appropriate RMS threshold value in a particular application was not described.

Recently, Maheshwari et al. described a multi-objective optimization (MOO) formulation for optimizing the design of the experiment using a combination of FIM-based metric and parameter correlation [15]. Because parameter correlations could not account for model nonlinearity, the strategy has the same drawback as FIM-based methods when applied to nonlinear models. In this work, we proposed a MOO MBDOE method using a combination of a FIM criterion and model curvature. We demonstrated the advantages of the proposed MOO MBDOE over FIM-based and other curvature-based methods in an application to the kinetic modeling of the fed-batch fermentation of baker’s yeast [30,31].

2. Model-Based Optimal Design of Experiments

We assume that the experimental data $\mathbf{y}\in {\mathrm{I}\mathrm{R}}^n$ are contaminated by additive random noise, as follows:

$$\mathbf{y}=\mathbf{\mu }+\mathbf{ϵ}$$

(1)

where $\mathbf{\mu }$ and $\mathbf{ϵ}$ denote the mean of the measurement data and the random noise, respectively. When the total number of data points n is greater than the number of parameters p, $\mathbf{\mu }$ spans a p-dimensional space $\mathsf{\Omega }\subset {\mathrm{I}\mathrm{R}}^n$, where

$$\mathsf{\Omega }=\{\mathbf{\mu }:\mathbf{\mu }=\mathbf{F}(\mathbf{x},\mathbf{u},\mathbf{\theta }),\mathbf{\theta }\in \mathsf{\Theta }\subset {\mathrm{I}\mathrm{R}}^p\}$$

(2)

Here, $\mathbf{x}\in {\mathrm{I}\mathrm{R}}^s$ denotes the state vector, $\mathbf{\theta }\in {\mathrm{I}\mathrm{R}}^p$ denotes the parameter vector, $\mathbf{u}\in {\mathrm{I}\mathrm{R}}^m$ denotes the input and $\mathbf{F}(\mathbf{x},\mathbf{u},\mathbf{\theta })$ denotes the vector of nonlinear model equations. The subspace $\mathsf{\Omega }$ is also called the expectation surface or the solution locus. For a dynamic system, the state $\mathbf{x}$ is often described by a set of ordinary differential equations (ODEs):

$$\frac{d\mathbf{x}(\mathbf{\theta },t)}{dt}=\mathbf{g}(\mathbf{x}(\mathbf{\theta },t),\mathbf{u},\mathbf{\theta }),\mathbf{x}(\mathbf{\theta },0)={\mathbf{x}}_0$$

(3)

The estimation of model parameters $\mathbf{\theta }$ from a given set of data $\mathbf{y}$ is typically formulated as a minimization of the weighted sum of squares of the difference between the model prediction $\mathbf{F}(\mathbf{x},\mathbf{u},\mathbf{\theta })$ and the measurement data $\mathbf{y}$. For example, the maximum likelihood estimator (MLE) of the model parameters for normally distributed data with known variance $\mathbf{V}$ is given by the minimum of the following objective function:

$$\mathsf{\Phi }\left(\mathbf{\theta }\right)={\left[\mathbf{y}-\mathbf{F}(\mathbf{x},\mathbf{u},\mathbf{\theta })\right]}^T{\mathbf{V}}^{-1}\left[\mathbf{y}-\mathbf{F}(\mathbf{x},\mathbf{u},\mathbf{\theta })\right]$$

(4)

When the model is a linear function of the parameters $\mathbf{F}(\mathbf{x},\mathbf{u},\mathbf{\theta })=\mathbf{X}\mathbf{\theta }$, $\mathbf{X}\in {\mathrm{I}\mathrm{R}}^{n\times p}$, then the parameter estimates are given by $\widehat{\mathbf{\theta }}={\left({\mathbf{X}}^T{\mathbf{V}}^{-1}\mathbf{X}\right)}^{-1}{\mathbf{X}}^T{\mathbf{V}}^{-1}\mathbf{y}$. In this case, the MLE is the minimum variance unbiased estimator of $\theta$, where the covariance matrix of the parameter estimates is given by ${\mathbf{V}}_{\mathbf{\theta }}={\left({\mathbf{X}}^T{\mathbf{V}}^{-1}\mathbf{X}\right)}^{-1}$. When the model is nonlinear (with respect to the parameters), the parameter estimates $\widehat{\mathbf{\theta }}=arg\min \mathsf{\Phi }\left(\mathbf{\theta }\right)$ do not necessarily correspond to the minimum variance estimator. According to the Cramér-Rao inequality [9], the inverse of the FIM provides a lower bound for the covariance of the parameter estimates $\widehat{\mathbf{\theta }}$, that is

$${\mathbf{V}}_{\mathbf{\theta }}\ge {\mathrm{FIM}}^{-1}={\left({\widehat{\dot{\mathbf{F}}}}^T{\mathbf{V}}^{-1}\widehat{\dot{\mathbf{F}}}\right)}^{-1}$$

(5)

where $\widehat{\dot{\mathbf{F}}}=\dot{\mathbf{F}}(\widehat{\mathbf{\theta }},\mathbf{x})=\frac{\partial \mathbf{F}(\mathbf{x},\mathbf{u},\mathbf{\theta })}{\partial \mathbf{\theta }}{|}_{\mathbf{\theta }=\widehat{\mathbf{\theta }}}$ is the first-order sensitivity matrix of $\mathbf{F}(\mathbf{x},\mathbf{u},\mathbf{\theta })$ with respect to the parameters $\theta$.

On the basis of the Cramér-Rao inequality, the FIM has been commonly used as a criterion of data informativeness in MBDOEs. Many methods for MBDOEs, such as those listed in

Table 1. Model-based designs of experiments (MBDOEs) using the Fisher information matrix (FIM).

FIM-Based MBDOE	Criterion
D-optimal	max ${\prod }_i{\lambda }_i$
A-optimal	max ${\sum }_i{\lambda }_i$
E-optimal	max $min\left({\lambda }_i\right)$
Modified E-optimal	max $\frac{min\left({\lambda }_i\right)}{max\left({\lambda }_i\right)}$

, are based on finding experimental conditions that optimize a FIM-based information metric. As shown in Equation (5), the FIM relies on a linearization of the model behavior with respect to the parameters. Essentially, the linearization replaces the expectation surface $\mathsf{\Omega }$ by its tangent plane at $\widehat{\mathbf{\theta }}$. The performance of the experimental design using a FIM-based criterion would therefore depend on whether (1) the model outputs vary proportionally with the parameter values (planar assumption), and (2) whether this proportionality is constant (uniform coordinate assumption) [32]. When the model is highly nonlinear with respect to the parameters, FIM-based MBDOEs may produce suboptimal designs [33,34]. A recent MOO MBDOE using a combination of a FIM criterion and parameter correlation has been shown to provide an improvement over FIM-based MBDOE methods [15]. However, this method also relies on the first-order parametric sensitivity matrix, and thus it could not account for model nonlinearity.

Curvature-based design of experiment methods such as the Q-optimality have been introduced to account for model nonlinearity by employing a second-order approximation of the model output. Here, the curvature of the expectation surface $\mathsf{\Omega }$ is captured using the second-order sensitivities of $\mathbf{F}(\mathbf{x},\mathbf{u},\mathbf{\theta })$ based on the Taylor series expansion:

$$\mathbf{F}(\mathbf{x},\mathbf{u},\mathbf{\theta })=\mathbf{F}(\mathbf{x},\mathbf{u},\widehat{\mathbf{\theta }})+\widehat{\dot{\mathbf{F}}}(\mathbf{\theta }-\widehat{\mathbf{\theta }})+\frac{1}{2}{(\mathbf{\theta }-\widehat{\mathbf{\theta }})}^T\widehat{\ddot{\mathbf{F}}}(\mathbf{\theta }-\widehat{\mathbf{\theta }})+\mathrm{O}\left({(\mathbf{\theta }-\widehat{\mathbf{\theta }})}^3\right)$$

(6)

where ${\widehat{\ddot{\mathbf{F}}}}_{ijk}=\frac{{\partial }^2{\mathrm{F}}_i(\mathbf{x},\mathbf{u},\mathbf{\theta })}{\partial {\theta }_j\partial {\theta }_k}{|}_{\mathbf{\theta }=\widehat{\mathbf{\theta }}}$ is the $n\times p\times p$ Hessian matrix. As mentioned in the Introduction, several curvature-based MBDOE methods are available, for example, by minimizing curvature or using a curvature threshold [30]. In this work, we employed a MOO approach based on curvatures for designing optimal experiments. The basic premise of our MBDOE is to select experimental conditions that maximize the informativeness of data and ensure that the model behaves relatively linearly with respect to the parameters. More specifically, our MBDOE uses two objective functions, the first of which involves the maximization of a FIM-based information metric, and the second of which involves the minimization of relative curvature measures [28]. The second objective function ensures that the FIM can provide a reliable measure of data informativeness.

2.1. Multi-Objective Design of Experiments Based on Curvatures

In this section, we derive the relative curvature measures by following the work of Bates and Watts [28]. We consider an arbitrary straight line in the parameter space passing through $\widehat{\mathbf{\theta }}$:

$$\mathbf{\theta }\left(b\right)=\widehat{\mathbf{\theta }}+b\mathbf{h}$$

(7)

where $\mathbf{h}=[h_1,h_2,\dots ,h_p]$ is a non-zero vector. As the scalar parameter b varies, a curve is traced through the expectation surface, also referred to as the lifted line, according to

$${\mathbf{\mu }}_{\mathbf{h}}\left(b\right)=\mathbf{\mu }(\widehat{\mathbf{\theta }}+b\mathbf{h})$$

(8)

The tangent line of this curve at $b=0$ is given by

$$\begin{array}{cc}\hfill {\mathbf{\mu }}_{\mathbf{h}}& ={\left[\frac{d{\mathbf{\mu }}_{\mathbf{h}}\left(b\right)}{db}\right]}_{\mathbf{\theta }=\widehat{\mathbf{\theta }},b=0}\hfill \\ & ={\left[\sum _{r=1}^p\frac{\partial \mathbf{F}(\mathbf{x},\mathbf{u},\mathbf{\theta })}{\partial {\theta }_r}\frac{\partial {\theta }_r\left(b\right)}{\partial b}\right]}_{\mathbf{\theta }=\widehat{\mathbf{\theta }},b=0}\hfill \\ & =\widehat{\dot{\mathbf{F}}}\mathbf{h}\hfill \end{array}$$

(9)

The set of all such tangent lines, that is, the column space of $\widehat{\dot{\mathbf{F}}}$, describes the tangent (hyper)plane at $\mathbf{\mu }\left(\widehat{\mathbf{\theta }}\right)$.

Meanwhile, the curvature measures come from a quadratic approximation of $\mathbf{\mu }$. In this case, the acceleration of $\mathbf{\mu }\left(b\right)$ at b = 0 can be written as follows:

$${\ddot{\mathbf{\mu }}}_{\mathbf{h}}={\mathbf{h}}^T\widehat{\ddot{\mathbf{F}}}\mathbf{h}=\sum _{i=1}^p\sum _{j=1}^p\frac{{\partial }^2\mathbf{F}(\mathbf{x},\mathbf{u},\mathbf{\theta })}{\partial {\theta }_i\partial {\theta }_j}{h}_i{h}_j$$

(10)

The acceleration vector ${\ddot{\mathbf{\mu }}}_{\mathbf{h}}$ can be subsequently decomposed into two components:

$${\ddot{\mathbf{\mu }}}_{\mathbf{h}}={\ddot{\mathbf{\mu }}}_{\mathbf{h}}^t+{\ddot{\mathbf{\mu }}}_{\mathbf{h}}^n$$

(11)

where at $\mathbf{\mu }\left(\widehat{\mathbf{\theta }}\right)$, ${\ddot{\mathbf{\mu }}}_{\mathbf{h}}^t$ is tangential to the tangent plane and ${\ddot{\mathbf{\mu }}}_{\mathbf{h}}^n$ is normal to the tangent plane. The tangential acceleration ${\ddot{\mathbf{\mu }}}_{\mathbf{h}}^t$ is also called the parameter-effect curvature [28] and provides a measure of nonlinearity along the parameter vector $\mathbf{h}$. The degree of the parameter-effect curvature can change upon reparameterization of the model. Meanwhile, the normal acceleration ${\ddot{\mathbf{\mu }}}_{\mathbf{h}}^n$ does not vary with model parameterization, and hence it is called the intrinsic curvature. Finally, the relative curvature measures in the direction of $\mathbf{h}$ are given by [28,32]:

$$K_{\mathbf{h}}^t=\frac{\parallel {\ddot{\mathbf{\mu }}}_{\mathbf{h}}^t\parallel }{\parallel {\ddot{\mathbf{\mu }}}_{\mathbf{h}}^t{\parallel }^2}$$

(12)

$$K_{\mathbf{h}}^n=\frac{\parallel {\ddot{\mathbf{\mu }}}_{\mathbf{h}}^n\parallel }{\parallel {\ddot{\mathbf{\mu }}}_{\mathbf{h}}^n{\parallel }^2}$$

(13)

Below, we describe the decomposition of the Hessian into the tangential and the normal component. We consider the QR-factorization of the Jacobian $\widehat{\dot{\mathbf{F}}}$, that is, $\widehat{\dot{\mathbf{F}}}=\mathbf{Q}\mathbf{R}=\mathbf{Q}\left[\begin{array}{c}\tilde{\mathbf{R}}\\ \mathbf{0}\end{array}\right]$. By rotating the parameter axes $(\mathbf{\theta }-\widehat{\mathbf{\theta }})$ into $\mathbf{\phi }=\tilde{\mathbf{R}}(\mathbf{\theta }-\widehat{\mathbf{\theta }})$, a new Jacobian matrix $\dot{\mathbf{U}}=\frac{d\mathbf{F}(\mathbf{x},\mathbf{u},\mathbf{\phi })}{d\mathbf{\phi }}{|}_{\mathbf{\phi }=0}$ can be computed as $\dot{\mathbf{U}}=\widehat{\dot{\mathbf{F}}}{\tilde{\mathbf{R}}}^{-1}$, which comprises the first p column vectors of $\mathbf{Q}$ (i.e., $\mathbf{Q}=\left[\begin{array}{cc}\dot{\mathbf{U}}& \mathbf{N}\end{array}\right]$). The remaining column vectors of $\mathbf{Q}$ (i.e., $\mathbf{N}$) are orthonormal to the tangent surface at $\mathbf{\phi }=0$. In the same manner, the Hessian matrix in the rotated axes can be written as $\ddot{\mathbf{U}}={\mathbf{L}}^T\widehat{\ddot{\mathbf{F}}}\mathbf{L}$, where $\mathbf{L}={\tilde{\mathbf{R}}}^{-1}$ and ${\ddot{U}}_{ijk}=\frac{{\partial }^2{F}_i(\mathbf{x},\mathbf{u},\mathbf{\phi })}{\partial {\phi }_j{\phi }_k}{|}_{\mathbf{\phi }=0}$. The decomposition of the Hessian into the tangential and normal components is given by the following equation [28]:

$$\ddot{\mathbf{A}}={\mathbf{Q}}^T\ddot{\mathbf{U}}={\left[\begin{array}{cc}\dot{\mathbf{U}}& \mathbf{N}\end{array}\right]}^T\ddot{\mathbf{U}}=\left[\begin{array}{cc}{\ddot{\mathbf{A}}}^t& {\ddot{\mathbf{A}}}^n\end{array}\right]$$

(14)

The matrices ${\ddot{\mathbf{A}}}^t$ and ${\ddot{\mathbf{A}}}^n$ respectively correspond to the parameter-effect and intrinsic curvature components of the Hessian.

To normalize the relative curvatures in Equations (12) and (13), Bates and Watts [28] used the scaling factor $\rho$, where $\rho =s\sqrt{p}$ and $s^2=\frac{{(\mathbf{y}-\widehat{\mathbf{\mu }})}^T(\mathbf{y}-\widehat{\mathbf{\mu }})}{n-p}$. Following the same procedure, we define the normalized relative curvatures as follows:

$${\gamma }_{\mathbf{h}}^t=\rho K_{\mathbf{h}}^t$$

(15)

$${\gamma }_{\mathbf{h}}^n=\rho K_{\mathbf{h}}^n$$

(16)

In addition, recasting $\mathbf{h}$ in the rotated axes as $\mathbf{h}=\mathbf{L}\mathbf{d}$, the tangent line ${\dot{\mathbf{\mu }}}_{\mathbf{L}\mathbf{d}}$ will have a unit norm (i.e., $\parallel {\dot{\mathbf{\mu }}}_{\mathbf{L}\mathbf{d}}\parallel =1$) when $\mathbf{d}$ is a unit vector. The computation of ${\gamma }_{\mathbf{h}}^t$ and ${\gamma }_{\mathbf{h}}^n$ is thus simplified into

$${\gamma }_{\mathbf{L}\mathbf{d}}^t=\rho \parallel {\mathbf{d}}^T{\ddot{\mathbf{A}}}^t\mathbf{d}\parallel ,\forall \mathbf{d}:\parallel \mathbf{d}\parallel =1$$

(17)

$${\gamma }_{\mathbf{L}\mathbf{d}}^n=\rho \parallel {\mathbf{d}}^T{\ddot{\mathbf{A}}}^n\mathbf{d}\parallel ,\forall \mathbf{d}:\parallel \mathbf{d}\parallel =1$$

(18)

In the proposed experimental design, the maximum of these curvature measures are used, where

$${\gamma }_{max}^t={\displaystyle \underset{\parallel \mathbf{d}\parallel =1}{\max }}{\gamma }_{\mathbf{L}\mathbf{d}}^t$$

(19)

$${\gamma }_{max}^n={\displaystyle \underset{\parallel \mathbf{d}\parallel =1}{\max }}{\gamma }_{\mathbf{L}\mathbf{d}}^n$$

(20)

As mentioned above, in formulating the MOO for the design of experiments, two design criteria have been taken into account. The first is that the experiment should be designed to maximize the informativeness of the data for parameter estimation. In this case, we employ an information metric based on the FIM. Meanwhile, the second design criterion in the MOO aims to minimize both the parameter-effect and intrinsic curvatures. The MOO formulation offers certain advantages, for example, that there is no need to prioritize any one of the criteria beforehand. Instead, we generate the Pareto set or Pareto frontier representing the set of solutions for which we cannot improve the value of one objective function without negatively affecting the other(s) [35].

Considering the kinetic ODE model given in Equation (3), our multi-objective formulation using the D-optimal criterion is given by

$$\begin{array}{c}\hfill {\displaystyle \underset{{\mathbf{x}}_{\mathbf{0}},{\mathbf{t}}_{\mathbf{sp}},\mathbf{u}\left(t\right)}{\max }}\prod _i{\lambda }_i\\ \hfill {\displaystyle \underset{{\mathbf{x}}_{\mathbf{0}},{\mathbf{t}}_{\mathbf{sp}},\mathbf{u}\left(t\right)}{\min }}{\gamma }_{max}^t+{\gamma }_{max}^n\end{array}$$

(21)

subject to

$$\begin{array}{c}\hfill \frac{d\mathbf{x}(\widehat{\mathbf{\theta }},t)}{dt}=\mathbf{g}(\mathbf{x}(\widehat{\mathbf{\theta }},t),\mathbf{u},\widehat{\mathbf{\theta }})\\ \hfill \mathbf{x}(\widehat{\mathbf{\theta }},0)={\mathbf{x}}_0\\ \hfill {\mathbf{x}}_0^L\le {\mathbf{x}}_0\le {\mathbf{x}}_0^U\\ \hfill {\mathbf{u}}_j^L\le {\mathbf{u}}_j\le {\mathbf{u}}_j^U\end{array}$$

(22)

where ${\lambda }_i$ is the ith eigenvalue of the FIM (Equation (5)). The first objective function can be substituted with other FIM-based metrics (see Table 1). The parameter vector $\widehat{\mathbf{\theta }}$ is either an initial guess of the parameter values or the parameter estimates from the current iteration of an iterative model identification procedure [6]. The decision variables may include the initial condition of the states ${\mathbf{x}}_0$, the sampling time points of measurements ${\mathbf{t}}_{sp}$, and the dynamic input $\mathbf{u}\left(t\right)$. In the case study below, we considered a control vector parametrization (CVP) of the input $u_i\left(t\right)$ as illustrated in Figure 2.

2.2. Numerical Implementation of the Curvature-Based MOO Design

As described in the previous section, the parameter-effect and intrinsic curvatures require the computation of the first- and second-order model sensitivities. For the ODE model in Equation (3), the first-order sensitivities can be calculated according to

$$\widehat{\dot{\mathbf{F}}}=\dot{\mathbf{F}}(\widehat{\mathbf{\theta }},\mathbf{x})=\frac{\partial \mathbf{F}\left(\mathbf{x}\right(t,\mathbf{u},\mathbf{\theta }\left)\right)}{\partial \mathbf{x}}\frac{\partial \mathbf{x}(t,\mathbf{u},\mathbf{\theta })}{\partial \mathbf{\theta }/\mathbf{\theta }}{|}_{\widehat{\mathbf{\theta }}}$$

(23)

The sensitivities in the above equation are normalized with respect to the parameter values. The last term on the right-hand side is the first-order sensitivities of the ODE model, which obey the following differential equation:

$$\frac{d}{dt}\frac{\partial \mathbf{x}}{\partial \mathbf{\theta }}=\frac{\partial \mathbf{g}}{\partial \mathbf{x}}\frac{\partial \mathbf{x}}{\partial \mathbf{\theta }}+\frac{\partial \mathbf{g}}{\partial \mathbf{\theta }},\frac{\partial \mathbf{x}}{\partial \mathbf{\theta }}{|}_{t=0}=0$$

(24)

Here, we have assumed that ${\mathbf{x}}_0$ is not part of the parameter estimation, but such an assumption can be easily relaxed. In the case study, the sensitivities $\frac{\partial \mathbf{x}}{\partial \mathbf{\theta }}$ were computed by solving the ODE in Equation (24) simultaneously to Equation (3), following a procedure known as the direct differential method [36]. Meanwhile, the Hessian matrix was approximated using a finite-difference method, as follows:

$${\widehat{\ddot{\mathbf{F}}}}_{ijk}=\left\{\begin{array}{c}\frac{F_i(\mathbf{\theta }+\Delta {\theta }_j{\mathbf{e}}_j)-2F_i\left(\mathbf{\theta }\right)+F_i(\mathbf{\theta }-\Delta {\theta }_j{\mathbf{e}}_j)}{\Delta {\theta }_j^2/{\theta }_j^2},\mathrm{for}j=k\hfill \\ \frac{F_i(\mathbf{\theta }+\Delta {\theta }_j{\mathbf{e}}_j+\Delta {\theta }_k{\mathbf{e}}_k)-F_i(\mathbf{\theta }+\Delta {\theta }_j{\mathbf{e}}_j-\Delta {\theta }_k{\mathbf{e}}_k)-F_i(\mathbf{\theta }-\Delta {\theta }_j{\mathbf{e}}_j+\Delta {\theta }_k{\mathbf{e}}_k)+F_i(\mathbf{\theta }-\Delta {\theta }_j{\mathbf{e}}_j-\Delta {\theta }_k{\mathbf{e}}_k)}{(\Delta {\theta }_j/{\theta }_j)(\Delta {\theta }_k/{\theta }_k)},\mathrm{for}j\ne k\hfill \end{array}\right.$$

(25)

where ${\mathbf{e}}_j$ is the jth elementary vector and uses 1% parameter perturbations (i.e., $\Delta {\theta }_j/{\theta }_j$ = 0.01). The second-order sensitivities above are also normalized with respect to the parameter values.

Meanwhile, the curvature measures ${\gamma }_{max}^t$ and ${\gamma }_{max}^n$ in Equations (19) and (20) were calculated from the Hessian matrix using the alternating least squares (ALS) method [37], an algorithm created to find the maximum singular value ${\sigma }_{max}$ of a three-dimensional matrix. Based on the definitions in Equations (19) and (20), the maximum curvature measures can be determined by computing the maximum singular values of the matrices $\rho {\ddot{\mathbf{A}}}^t$ and $\rho {\ddot{\mathbf{A}}}^n$, respectively. More specifically, we implemented the ALS method to solve for

$${\sigma }_{max}\left(\mathbf{B}\right)=\underset{\parallel \mathbf{r}\parallel =\parallel \mathbf{s}\parallel =1}{\max }\sum _{i=1}^m{\mathbf{r}}_{\mathbf{i}}{\mathbf{s}}^T{\mathbf{B}}_{\mathbf{i}}\mathbf{s}$$

(26)

where $\mathbf{B}$ is either $\rho {\ddot{\mathbf{A}}}^t$ or $\rho {\ddot{\mathbf{A}}}^n$. The ALS algorithm started with initial guess values of the vectors $\mathbf{r}$ and $\mathbf{s}$ and used the above equation to solve for one variable while fixing the other in an alternating manner. Zhang and Golub showed that the method linearly converges in a neighbourhood of the optimal solution [37].

In the case study, the MOO problem was solved using the non-dominated sorting genetic algorithm II (NSGAII) in MATLAB, producing a Pareto frontier in the space of the objective functions [38]. We employed a population size of 300 and set the number of generations to 50 times the number of parameters (i.e., 1450). We recasted a maximization of an objective function as the minimization of its negative counterpart. The optimal design was selected from the Pareto frontier by balancing the trade-offs among the objective functions. More specifically, we first normalized the objective functions such that their values on the Pareto frontier ranged between 0 and 1. Finally, we chose among the solutions on the Pareto frontier that which minimized the Euclidean distance of all (normalized) objective functions as the final design.

3. Results

3.1. MBDOEs of Baker Yeast Fermentation Model

We evaluated the performance of the proposed MBDOE in an application to a kinetic model of a fed-batch fermentation of baker’s yeast [30,31]. In addition to a D-optimal criterion, we also implemented A-optimal, E-optimal and modified E-optimal criteria (see Table 1) with our MOO MBDOE. We compared the performance of our method to other MBDOEs, including (a) FIM-based MBDOEs, that is, D-optimal, A-optimal, E-optimal and modified E-optimal designs; (b) a D-optimal design with a curvature threshold [30]; (c) a Q-optimal MBDOE [29]; and (d) a MOO MBDOE using parameter correlation [15]. In total, we applied and compared 14 MBDOE methods. For the optimizations in (a), (b) and (c), we employed the enhanced scatter search metaheuristic (eSSm) algorithm [39,40,41]. For the MOO in (d), we used the optimization algorithm and optimal Pareto point selection, as described in the previous section.

In the fed-batch fermenter model, cellular growth and product formation are captured by the biomass variable $x_1$, which is assumed to rely on a single substrate variable $x_2$. The fermenter operates at a constant temperature and the feed is free from product. The model equations are given by

$$\begin{array}{c}\hfill \frac{d{x}_1}{dt}=(r-u_1-{\theta }_4)x_1,x_1\left(0\right)=x_{10}\\ \hfill \frac{d{x}_2}{dt}=-\frac{r{x}_1}{{\theta }_3}+u_1(u_2-x_2),x_2\left(0\right)=0.1\\ \hfill r=\frac{{\theta }_1{x}_2}{{\theta }_2+x_2}\end{array}$$

(27)

where the input $u_1$ is the dilution factor (in the range of 0.05–0.20 h${}^{-1}$) and the input $u_2$ is the substrate concentration in the feed (in the range of 5–35 g/L). In the model, the biomass growth follows Monod-type kinetics. The parameters ${\theta }_1$ and ${\theta }_2$ are the Monod kinetic parameters, ${\theta }_3$ is the yield coefficient, and ${\theta }_4$ is the cell death rate constant.

In the MBDOE, the design variables consisted of the initial condition of the biomass $x_1\left(0\right)$ in the range between 1 and 10 g/L, 10 measurement sampling times ($t_{sp}$), and the inputs $u_1\left(t\right)$ and $u_2\left(t\right)$. The piecewise-constant dynamic inputs were each parametrized using the CVP, as shown in Figure 2. Thus, the MOO was performed with 29 design parameters ($x_1\left(0\right)$, 10 $t_{sp}$’s, 10 $u_{i,j}$’s, and 8 $t_{sw}$’s). The length of the time interval between two successive measurement sampling points was constrained to be between 1 and 20 h, while that between two input switching times was bounded between 2 and 20 h. The calculations of the Jacobian and Hessian matrices in MBDOEs were made using parameter values ${\mathbf{\theta }}_d=[{\theta }_1,{\theta }_2,{\theta }_3,{\theta }_4]=[0.5,0.5,0.5,0.5]$ [15,42], which were different from the “true” parameter values used for noisy data generation in the next section. The reason for using a different parameter set in the MBDOE to the true values was to emulate the typical scenario in practice, for which one would start only with an estimate or guess of the model parameters. Figure 3 and Figure 4 show the optimal dynamic inputs and data sampling times resulting from all the MBDOE methods mentioned above (see also the Pareto frontiers in Figures S1 and S2 in the Supplementary Materials). Meanwhile,

Table 2. Optimal initial condition of biomass $x_1\left(0\right)$ (g/L) from model-based design of experiments (MOO: multi-objective optimization).

Design Criterion	${\mathit{x}}_1\left(0\right)$
D-optimal	10.0
MOO D-optimal and curvatures	10.0
MOO D-optimal and correlation	10.0
A-optimal	10.0
MOO A-optimal and curvatures	9.9
MOO A-optimal and correlation	10.0
E-optimal	10.0
MOO E-optimal and curvatures	10.0
MOO E-optimal and correlation	10.0
Modified E-optimal	10.0
MOO modified E-optimal and curvatures	10.0
MOO modified E-optimal and correlation	10.0
Threshold curvature	8.2
Q-optimal	5.5

gives the optimal initial biomass concentration $x_1\left(0\right)$.

3.2. Performance Evaluation

For each of the optimal experimental designs above, we generated in silico datasets by simulating the ODE model using the parameter values ${\mathbf{\theta }}^{\ast }$= [0.31, 0.18, 0.55, 0.05], as reported in previous publications [15,42]. We subsequently added independent and identically distributed (i.i.d.) Gaussian random white noise to the model simulations using a relative variance of 0.04 for both $x_1\left(t\right)$ and $x_2\left(t\right)$ [15,42]. For each in silico dataset, we then performed a parameter estimation using the resulting data ($y_1$ and $y_2$) by maximum likelihood estimation, that is, by minimizing

$$\mathsf{\Phi }\left(\theta \right)=\frac{1}{{\sigma }^2}\sum _{i=1}^{10}{[y_1\left(t_i\right)-x_1(t_i,\mathbf{\theta })]}^2+{[y_2\left(t_i\right)-x_2(t_i,\mathbf{\theta })]}^2$$

(28)

We also employed the following constraints for $\mathbf{\theta }$ in the optimization above:

$$0.05\le {\theta }_1,{\theta }_2,{\theta }_3\le 0.98$$

and

$$0.01\le {\theta }_4\le 0.98$$

Finding the globally optimal solution to the parameter estimation in Equation (4) is challenging. Here, we solved the constrained parameter optimization problem using the interior-point algorithm (implemented by the subroutine fmincon function in MATLAB) with the true parameter values ${\mathbf{\theta }}^{\ast }$ as the initial guess. By employing the true values as the initial starting point of the optimization, we expected that the parameter accuracy would mainly be affected by the experimental design and not by the ability of the parameter optimization algorithm to find the globally optimal solution.

We repeated the in silico data generation and parameter estimation as described above 100 times, which resulted in a set of 100 parameter estimates. The performance of each MBDOE was assessed by the average accuracy of the parameter estimates, measured by the average of the normalized mean-square error (nMSE):

$$\overline{\mathrm{nMSE}}=\frac{1}{4}\sum _{i=1}^4{\mathrm{nMSE}}_{\mathrm{i}}$$

(29)

where

$${\mathrm{nMSE}}_{\mathrm{i}}=\frac{\mathrm{variance}\left({\widehat{\theta }}_i\right)-{\mathrm{bias}}^2\left({\widehat{\theta }}_i\right)}{{\left({\theta }_i^{\ast }\right)}^2},i=1,2,3,4$$

(30)

The variance of ${\widehat{\theta }}_i$ was computed using the set of 100 parameter estimates, while the bias was calculated as the difference between the average of ${\widehat{\theta }}_i$ and ${\theta }_i^{\ast }$.

Table 3. Model-based design of experiment (MBDOE) performance on the fed-batch fermentation of baker’s yeast model. The overall parameter accuracy is represented by the average of the normalized mean-square error (nMSE). The reported parameter values and errors are the averages and standard deviations from 100 repeated runs of parameter estimation.

Design Criterion	$\overline{\mathrm{nMSE}}$	${\mathit{\theta }}_1$ ± ${\mathit{SD}}_{{\mathit{\theta }}_1}$	${\mathit{\theta }}_2$ ± ${\mathit{SD}}_{{\mathit{\theta }}_2}$	${\mathit{\theta }}_3$ ± ${\mathit{SD}}_{{\mathit{\theta }}_3}$	${\mathit{\theta }}_4$ ± ${\mathit{SD}}_{{\mathit{\theta }}_4}$
D-optimal	7.06 × 10${}^{-3}$	0.3107 ± 0.0102	0.1831 ± 0.0276	0.5505 ± 0.0125	0.0502 ± 0.0026
MOO D-optimal and curvatures	4.71 × 10${}^{-3}$	0.3099 ± 0.0056	0.1825 ± 0.0233	0.5496 ± 0.0099	0.0499 ± 0.0018
MOO D-optimal and correlation	5.36 × 10${}^{-3}$	0.3117 ± 0.0134	0.1781 ± 0.0151	0.5543 ± 0.0270	0.0508 ± 0.0049
A-optimal	2.35 × 10${}^{-1}$	0.3294 ± 0.0659	0.2399 ± 0.1387	0.5841 ± 0.1083	0.0558 ± 0.0181
MOO A-optimal and curvatures	1.42	0.3669 ± 0.0947	0.5267 ± 0.2230	0.5548 ± 0.1333	0.0510 ± 0.0244
MOO A-optimal and correlation	4.82	0.0863 ± 0.0499	0.8927 ± 0.2555	0.2879 ± 0.1928	0.0177 ± 0.0263
E-optimal	8.01 × 10${}^{-2}$	0.3180 ± 0.0420	0.2026 ± 0.0956	0.5473 ± 0.0159	0.0496 ± 0.0026
MOO E-optimal and curvatures	3.33 × 10${}^{-3}$	0.3083 ± 0.0095	0.1829 ± 0.0164	0.5502 ± 0.0183	0.0500 ± 0.0026
MOO E-optimal and correlation	8.19 × 10${}^{-3}$	0.3108 ± 0.0164	0.1824 ± 0.0213	0.5552 ± 0.0304	0.0509 ± 0.0055
Modified E-optimal	6.99 × 10${}^{-2}$	0.3137 ± 0.0165	0.1986 ± 0.0920	0.5498 ± 0.0144	0.0502 ± 0.0033
MOO modified E-optimal and curvatures	3.44 × 10${}^{-4}$	0.3095 ± 0.0036	0.1789 ± 0.0034	0.5491 ± 0.0073	0.0500 ± 0.0013
MOO modified E-optimal and correlation	2.27 × 10${}^{-3}$	0.3088 ± 0.0048	0.1820 ± 0.0160	0.5486 ± 0.0047	0.0496 ± 0.0013
Threshold curvature	1.29 × 10${}^{-2}$	0.3144 ± 0.0307	0.1857 ± 0.0339	0.5500 ± 0.0155	0.0502 ± 0.0032
Q-optimal	1.91 × 10${}^{-2}$	0.3085 ± 0.0178	0.1757 ± 0.0216	0.5514 ± 0.0236	0.0504 ± 0.0119

gives the average nMSE of the parameter estimates from each MBDOE under consideration.

4. Discussion

As shown in Figure 3, the MBDOEs prescribed manipulating the input $u_1\left(t\right)$ mostly at the beginning of the experiment and the input $u_2\left(t\right)$ for the entire duration of the experiment. For the majority of the MBDOEs in this study, the optimal sampling times spread unevenly over the duration of the experiment (see Figure 4). A more detailed comparison between Figure 3 and Figure 4 showed that the optimal sampling points were typically placed before and after a change in the dynamic inputs $u_1\left(t\right)$ and $u_2\left(t\right)$. The exception to this observation was for the optimal design using the A-optimal criterion, which gave the worst parameter accuracy among the MBDOEs considered.

The consideration of model curvature using the proposed MOO MBDOE generally led to improved parameter accuracy over using only model curvature (i.e., Q-optimal and threshold curvature) or using only FIM-based criteria. The lowest nMSE came from the MOO MBDOE design using the modified E-optimal with model curvature. In comparison to MOO MBDOE using parameter correlation, employing model curvature in the MOO framework gave better experimental designs with lower average nMSEs. Meanwhile, Q-optimality and curvature thresholding strategies provided better nMSEs than the majority of the FIM-based criteria, except the D-optimal design. Finally, the optimal experiments based on the A-optimal criterion, either alone or in MOO MBDOE, performed poorly. The poor performance of the A-optimal design has also been reported in a previous publication [15].

The obvious drawback of curvature-based MBDOEs in comparison to FIM-based strategies is the higher computational cost associated with computing the Hessian matrix. While the number of first-order sensitivities (Jacobian) increases linearly with the number of parameters p, the number of second-order sensitivities scales with $p^2$. Fortunately, the calculation of the Hessian matrix can be easily parallelized and implemented using multiple computing cores. In practice, one often focuses on only a subset of the model parameters, and therefore the MBDOE is typically done for a handful of parameters.

We note that the MBDOE methods considered in this work consider only parametric uncertainty in the models and assume that uncertainty in model equations, that is, structural uncertainty, is not significant. For certain types of models, such as generalized mass action and S-system models [43,44], model structural uncertainty can be treated as parametric uncertainty, and therefore the MBDOE strategies developed here could be applied. As mentioned in the Introduction, MBDOE methods for discriminating between model structures have been developed, many of which are based on the Bayesian approach. Furthermore, in applications for which there exists intrinsic parametric variability, for example, batch-to-batch variability in cell culture fermentation processes, Bayesian MBDOE methods would be more suitable than FIM-based strategies, as Bayesian methods are able to incorporate the prior (intrinsic) distribution of the parameter values in the design. Nevertheless, as demonstrated in the case study, even when the MOO MBDOEs were performed using model parameters that were quite different from the true values, the resulting optimal designs led to precise and accurate parameter estimates. Meanwhile, biological systems, like other complex systems, have been argued to be sloppy. In the context of our work, sloppy systems lead to mathematical models whose FIMs have eigenvalues that are logarithmically spread evenly over large orders of magnitude [45]. In other words, the system behavior is sensitive to or is controlled by a small number of parameter combinations (along the FIM eigenvectors corresponding to large eigenvalues). At the same time, there exist many parameter combinations that can be varied without affecting the system behavior. Such sloppiness could arise in a system governed by processes that span large and evenly distributed length and/or time scales, such that there exists no clear separation between relevant and irrelevant mechanisms. A recent study demonstrated that in the case of sloppy systems, reducing the model parametric uncertainty by MBDOEs beyond a certain point might not necessarily translate to any improvement in model prediction accuracy [45]. However, it is possible to construct reduced-order models of sloppy systems, whose parameters correspond to the important parameter combinations [46,47]. Parameter estimation and MBDOE strategies can then be applied to these reduced models.

5. Conclusions

Existing MBDOE methods for parameter estimation mostly rely on the FIM to define information criteria. Because the FIM is based on first-order sensitivities with respect to the model parameters, the related MBDOEs may perform poorly for nonlinear models. Here, a new MBDOE using a MOO framework was presented, employing the maximization of a FIM-based information metric and the minimization of model curvatures. The application to a model of the fermentation of baker’s yeast demonstrated that accounting model nonlinearity through model curvatures in designing the experiment could lead to improved parameter accuracy over using only a FIM-based criterion. The proposed MOO MBDOE also outperformed other curvature-based designs, including the Q-optimality and curvature thresholding and another MOO MBDOE strategy using parameter correlation. The use of the MOO framework further gives flexibility to accommodate other criteria that may arise in a particular application, in the design of experiments.