Optimal Experimental Design for Parameter Estimation of an IL-6 Signaling Model

IL-6 signaling plays an important role in inflammatory processes in the body. While a number of models for IL-6 signaling are available, the parameters associated with these models vary from case to case as they are non-trivial to determine. In this study, optimal experimental design is utilized to reduce the parameter uncertainty of an IL-6 signaling model consisting of ordinary differential equations, thereby increasing the accuracy of the estimated parameter values and, potentially, the model itself. The D-optimality criterion, operating on the Fisher information matrix and, separately, on a sensitivity matrix computed from the Morris method, was used as the objective function for the optimal experimental design problem. Optimal input functions for model parameter estimation were identified by solving the optimal experimental design problem, and the resulting input functions were shown to significantly decrease parameter uncertainty in simulated experiments. Interestingly, the determined optimal input functions took on the shape of PRBS signals even though there were no restrictions on their nature. Future work should corroborate these findings by applying the determined optimal experimental design on a real experiment.


Introduction
Mathematical models of intracellular signaling pathways are important for understanding and predicting how cells respond to certain stimuli.Such models can be modified readily when new findings become available, and can be a useful tool in directing new studies based on hypotheses generated by the model's predictions.
Interleukin-6 (IL-6) is a cytokine involved in a variety of inflammatory processes [1][2][3][4].Understanding the signaling pathways associated with extracellular IL-6 excitation is important for elucidating and modulating the biological response to inflammation.Because chronic inflammation can cause tissue damage and poses a serious health risk during chronic infection or autoimmune conditions [5][6][7], therapeutic treatments for chronic inflammation are an active area of research.These efforts can be augmented through the use of mathematical models for inflammatory signaling, where IL-6 plays a major role.
An IL-6 signaling model (Appendix A) was previously developed in [1][2][3][4].Originally, a detailed model containing 77 ordinary differential equations (ODEs) and 128 parameters was derived, followed by a model simplification procedure to decrease the number of ODEs and parameters to 13 and 19, respectively.The simplified model contains only the variables and parameters deemed necessary for representing the dynamics of the signaling system, as determined by parameter sensitivity analysis and observability analysis [2].This simplified model is the subject of the present study, and will henceforth be referred to as the 'IL-6 model' or the 'model' while the original detailed model will be referred to as the 'original IL-6 model' or the 'original model'.
In this study, optimal experimental design was applied to increase the accuracy of the 19 parameters of the IL-6 model.The parameters in this model correspond to rate constants for the chemical reactions in the signaling pathway represented by the model.By increasing the accuracy-and conversely reducing the uncertainty-of the model parameters, it is possible to obtain a more accurate model through iterative model adjustment.To demonstrate and examine this step in the model development process, optimal experimental design methods were applied in this study to minimize the uncertainty of the IL-6 model parameters estimated from least squares with experimental data.
Optimal experimental design has been utilized for decades in a variety of settings in which it is of interest to maximize efficiency of resource use and obtain a significant amount of information from experiments with acceptable cost [8][9][10][11][12][13][14][15][16].Recently, as biological modeling and systems biology have emerged as an important area in biomedical research, optimal experimental design applied to biological experimental systems has become more popular [17][18][19][20][21][22][23][24][25][26][27][28]; additionally, optimal experimental design has been recognized as a valuable tool in optimal control for several decades [29].For example, Jones et al. [13] maximized production of an exogenous commodity chemical in metabolically engineered E. coli using an empirical modeling method similar to those used in [15,16] to maximize the efficacy of drug delivery.Weber [26] utilized optimal experimental design to maximize model prediction accuracy for a model of vesicle transport via the trans-Golgi network.Bandara [18] performed optimal experimental design to reduce parameter uncertainty in a model of phosphatidylinositol 3,4,5-trisphosphate signaling.These studies demonstrate the effectiveness of optimal experimental design for obtaining maximally informative experimental data.
Here, optimal experimental design is applied to the problem of maximizing parameter accuracy for the IL-6 model.In particular, the D-optimality criterion, applied to the Fisher information matrix (FIM), is maximized over a set of IL-6 concentration input functions to determine an optimal dynamic IL-6 input profile for exciting the signaling system to generate data for least squares estimation of the model parameters [30].The experimental design constraints considered in the problem are based on available resources and limitations present in a typical laboratory capable of performing the designed experiments.Namely, it is assumed that (a) the measurements of protein concentrations for two transcription factors are recorded as a time series; (b) the sampling time between measurements does not change during the course of an experiment; (c) the model represents signal transduction in rat hepatocytes which are stimulated in vitro with IL-6 following a dynamic input function; and (d) the IL-6 concentration is kept below cytotoxic levels.A piecewise constant IL-6 input function was computed by solving the optimal experimental design problem with the D-optimality criterion, operating on the Fisher information matrix, as the objective function.Since the Fisher information matrix contains only local information, a sensitivity matrix was also computed using the Morris method and the D-optimality criterion was applied to the sensitivity matrix as well in order to corroborate that the results are not just local in nature.The optimal IL-6 input function was found to substantially decrease the model parameter uncertainty in simulated least squares fits compared to a constant stimulation with IL-6.
The paper is organized as follows: optimal experimental design and its application to ODE models are presented in Section 2; application of optimal experimental design to the specific case of the IL-6 model is presented in Section 3 and results from solving the optimal experimental design problem are presented in Section 4, including simulation results for calculating parameter uncertainty; Section 5 discusses implications and suggests further applications of the optimal experimental design methodology presented here.The IL-6 model is included for reference in Appendix A.

Optimal Experimental Design
An ordinary differential equation model, e.g., for describing signal transduction, can be written as where x is a time-dependent vector of state variables, u is a time-dependent controlled input to the system, p is a vector of constant parameters, and y is a vector of measured quantities related to the model variables.Often, as is the case here, y is simply a subset of state variables from the vector x that are measured experimentally over time.While these measurements y are discrete in nature, the assumption that y is a continuous variable can be made if the sampling frequency is sufficiently high.Optimal experimental design involves maximizing a criterion function that indicates the quantity of information gained by a given experiment, often in the context of model identification [8][9][10][11][12]14,[17][18][19][20][21][22][23]29,31].Several commonly used criterion functions for experimental design exist.One of these is A-optimality, which seeks to minimize the trace of the inverse of the Fisher information matrix.This criterion results in minimizing the average variance of the estimates of the regression coefficients.Another popular approach is E-optimality which maximizes the smallest eigenvalue of the Fisher information matrix.However, the most popular approach is D-optimal experimental design, i.e., which maximizes the determinant of the FIM [14,19,24].The D-optimality criterion was chosen for this work as it seeks to minimize the covariance of the parameter estimates for a specified model [14].It should be noted that no experimental design is optimal in all aspects and maximizing one particular criterion often negatively affects the experimental design of other criteria.That being said, the D-optimality criterion has found widespread use in practice as it results in experimental designs that have many of the properties that one usually looks for.For numerical considerations, the natural logarithm of the determinant of the FIM was taken for the problem presented.The optimal experimental design objective can be written generally as where F is the FIM, defined in more detail below, and ϕ D is the D-optimality criterion.
Taking the determinant of the parameter covariance matrix corresponds to computing the volume of the parameter space that would allow for a solution to the least squares parameter fitting problem [14,18].The volume of this parameter space represents the covariance, or uncertainty, of the parameters.Because the FIM is related to the inverse of the covariance matrix [8,10], maximizing the determinant of the FIM results in minimizing the determinant of the covariance matrix, and thus minimizing the volume of the parameter space, or minimizing parameter uncertainty.
The FIM is computed from the sensitivity coefficients of the model, which can be an ODE model as shown in Equation (1).Local sensitivity coefficients are defined as where x i is the i-th state variable in the model, p j is the j-th model parameter in the vector p, and t is the time at which the partial derivative is evaluated [2,3,32].The local sensitivity coefficient represents the change in a model state with respect to a change in the value of a parameter, and is a function of time and the parameter vector.At every time point during an experiment, or a model simulation, a sensitivity coefficient can be calculated for any of the state variables.However, the FIM contains the sensitivity coefficients for only those state variables which are measured in experiments for generating data to be used in parameter estimation using least squares fitting, i.e., the sensitivity coefficients for the vector y.The FIM is written as where S is the sensitivity matrix [32] and n p is the number of parameters in the model.Here, the formula shown for S corresponds to an experiment in which two state variables are measured in a time series from time t 1 to t N ; one can extend the formula to the general case of any number of experimentally measured outputs, however, two outputs are realistic for the application investigated in this paper.The local sensitivity coefficients for all of the state variables of the model can be calculated by solving the system of ordinary differential equations d dt where x is the column vector of state variables, p is the column vector of model parameters, and f is the column vector of functions defining the model, as in Equation ( 1) [32].To calculate the D-optimality criterion value for a given experiment design, the model ODEs (Equation ( 1)) are solved simultaneously with the sensitivity equations (Equation ( 6)), the FIM is constructed (Equations ( 4) and ( 5)), and the determinant can then be computed.To determine an optimal experiment, an optimization problem is solved which searches through different experimental designs, computing the D-optimality criterion in this way for every iteration.The Morris method [3] can also be utilized for calculating sensitivity coefficients in a covariance matrix.It should be noted that such an approach results in more than local information due to the certain properties of the Morris method.Using the Morris method, global sensitivity coefficients are calculated by repeatedly sampling parameter values from a distribution.Finite difference approximations of the local sensitivities are then computed at each value of the sampled parameter vector, and the finite difference approximations are averaged to obtain each global sensitivity coefficient as shown in Equations ( 7) and (8).
Here, d ijk is a finite difference derivative approximation for the local sensitivity of y i with respect to parameter p j at the k-th sampled value of p j , ∆ jk defines the sampled value of p j depending on the distribution of p j , and n d is the number of samples chosen for p j .This method for calculating sensitivity coefficients takes into account the uncertainty of the parameter values by averaging over samples from the parameter distribution.Therefore, using the Morris method for calculating the sensitivity coefficients in the FIM does not rely on knowledge of an exact value for the parameters.In fact, the exact values of the parameters are by definition unknown-otherwise, one would not need to estimate the parameters.For this problem, the parameter distributions were considered to be normal with a mean at the nominal parameter values and a standard deviation of one-tenth the value of each parameter.

Formulation of the Optimal Experimental Design Problem for the IL-6 Model
The values of the parameters p in ODE models are generally unknown and must be estimated using experimentally measured data points [20,30,33].The problem considered here is to optimize the function for the controlled variable u in order to minimize the uncertainty of the parameters estimated by least squares parameter estimation.
In the IL-6 model, u represents the IL-6 concentration in the local environment of the cells being stimulated by IL-6.The IL-6 concentration at a given time determines the signaling behavior of the cells, which is characterized by the concentrations of each signaling molecule in the pathway at time t.Thus by modulating the input IL-6 concentration, u, as a function of time, the signaling dynamics can be influenced so that measurements of components of the signaling pathway can yield maximal information about the model parameters.
The measured state variables for this problem were chosen to be STAT3N*-STAT3N* and C/EBPβ, two transcription factors of the IL-6 signaling pathway which are directly involved in transcribing DNA to RNA.These measurements were chosen because of the availability of a fluorescent reporter system for measuring the concentrations of these two proteins [2,34].To design an optimal experiment for minimizing parameter uncertainty, experiments were considered in which the concentrations of these two proteins are measured in a time series every 45 min for 22 h using a fluorescent reporter system.This time sequence for image acquisition was utilized previously to obtain initial estimates for the parameters of the model [2].However, for the initial parameter estimation, a constant IL-6 concentration at an arbitrary level of 100 ng/mL was utilized as the input function [2].In the present work, optimal experimental design was applied to optimize the input function over a continuous set of time-dependent input functions while utilizing the same time sequence for data acquisition as was utilized in [2].This allows for evaluating whether the parameter uncertainty can be decreased by optimizing the input function while holding other experimental control decisions constant.
For this problem, the set of input functions to optimize was chosen to be the piecewise constant input functions with fixed and equal time intervals and IL-6 concentrations bounded between 0 and 7.5 nM.These input functions were chosen because they can be implemented experimentally due to the long intervals during which a concentration is constant, and because they can be parameterized so as to arrive at an optimization problem with a finite number of variables.Specifically, the number of optimization variables is the number of concentration levels allowed in the piecewise constant function.This is illustrated in Figure 1 with a piecewise constant input function for IL-6 that has r = 4 concentration levels.As an example, to determine the optimal input function with r = 4 concentration levels, the optimal experimental design problem would be solved with four optimization variables representing the four concentration levels of such an input function.
These piecewise constant input functions can be written as where c k is the k-th concentration level in the vector c, 'step' is the Heaviside step function, r is the number of concentration levels in the input function, and ∆t = 22 h/r is the time interval for each concentration level in the input function.
By changing the vector c to modulate the IL-6 input function, the solution of the model ODEs and local sensitivity ODEs are modulated.This causes the FIM, and thus the D-optimality criterion value, to be a function of c.Therefore, the optimal experimental design problem for minimizing parameter uncertainty in the IL-6 model can be written as max 0≤c≤7.5 nM ln|F(c)|, (10) where |•| is the determinant and F is calculated from Equations ( 4) and ( 5) with y 1 from Equation (5)  being the concentration of STAT3N*-STAT3N*, y 2 being the concentration of C/EBPβ, and p being the parameters of the IL-6 model.Specifically, in order to calculate the FIM, the IL-6 signaling ODE model is numerically integrated simultaneously with the local sensitivity equations (Equation ( 6)) for a given input function represented by c, and the FIM is constructed according to Equations ( 4) and ( 5).This ODE integration is carried out in every iteration of the optimization problem to evaluate the D-optimality criterion value until the optimization solver determines a solution.
Processes 2017, 5, 49 6 of 12 ODE model is numerically integrated simultaneously with the local sensitivity equations (Equation ( 6)) for a given input function represented by c, and the FIM is constructed according to Equations ( 4) and ( 5).This ODE integration is carried out in every iteration of the optimization problem to evaluate the D-optimality criterion value until the optimization solver determines a solution.When solving the problem using the Morris method to obtain global sensitivity coefficients, the local sensitivity coefficients in the FIM are replaced by the global sensitivity coefficients calculated by the Morris method.For the Morris method, the local sensitivity ODEs do not need to be solved; rather, the model ODEs are solved repeatedly using the sampled parameter values, and the finite difference approximations of the partial derivatives are computed and averaged.
To solve the optimal experimental design problem for each value of r from 1 to 6, the MATLAB function 'fminsearch' was utilized [35].Additionally, to account for the bound constraints on the concentration levels, a stop flag was imposed in the program for running the optimization solver.Multiple initial guesses were utilized for each value of r to avoid local optima (between 5 and 18 initial guesses were used depending on r, in order to cover the range of possible qualitative input function shapes for each value of r).

Optimal Experimental Design for the IL-6 Signaling Model
Solving the optimal experimental design problem for minimizing parameter uncertainty in the IL-6 model for values of r from 1 to 6 resulted in an optimal IL-6 input function for each value of r.The optimal solutions for each r are listed in Table 1, which also lists several sub-optimal input functions for comparison.Note that the D-optimality criterion value is greatest for the optimal input function with r = 6 (see columns 1-4 in Table 1).The range of D-optimality criterion values for the input functions listed in the table (column 4) is very wide if one considers that the values given are the logarithm of the determinant of the FIM rather than the determinant itself.As expected, the D-optimality criterion values for the optimal input functions increase with r (see Table 1, column 4); this is due to the degrees of freedom added when r is increased.When solving the problem using the Morris method to obtain global sensitivity coefficients, the local sensitivity coefficients in the FIM are replaced by the global sensitivity coefficients calculated by the Morris method.For the Morris method, the local sensitivity ODEs do not need to be solved; rather, the model ODEs are solved repeatedly using the sampled parameter values, and the finite difference approximations of the partial derivatives are computed and averaged.
To solve the optimal experimental design problem for each value of r from 1 to 6, the MATLAB function 'fminsearch' was utilized [35].Additionally, to account for the bound constraints on the concentration levels, a stop flag was imposed in the program for running the optimization solver.Multiple initial guesses were utilized for each value of r to avoid local optima (between 5 and 18 initial guesses were used depending on r, in order to cover the range of possible qualitative input function shapes for each value of r).

Optimal Experimental Design for the IL-6 Signaling Model
Solving the optimal experimental design problem for minimizing parameter uncertainty in the IL-6 model for values of r from 1 to 6 resulted in an optimal IL-6 input function for each value of r.
The optimal solutions for each r are listed in Table 1, which also lists several sub-optimal input functions for comparison.Note that the D-optimality criterion value is greatest for the optimal input function with r = 6 (see columns 1-4 in Table 1).The range of D-optimality criterion values for the input functions listed in the table (column 4) is very wide if one considers that the values given are the logarithm of the determinant of the FIM rather than the determinant itself.As expected, the D-optimality criterion values for the optimal input functions increase with r (see Table 1, column 4); this is due to the degrees of freedom added when r is increased.To test whether higher D-optimality criterion values for an input function corresponded to lower uncertainty, simulations were run for each of the input functions listed in Table 1.For these simulations, data were generated by utilizing the original IL-6 model and adding normally distributed noise to the measurements.Specifically, for a given input function, a simulation was run by integrating the original IL-6 model, adding Gaussian noise with a mean of 0 and a standard deviation of 1 at the measurement time points (every 45 min for 22 h), and fitting the parameters of the simplified IL-6 model to the simulated data.The standard deviation of 1 was chosen to provide a reasonable signal to noise ratio for the model variables.For each input function, 30 simulations were run and parameters were fitted for each simulation.A parameter covariance matrix was then calculated for each of the input functions using the parameter fits from the 30 simulations that were run for each input function.The norm and trace of these covariance matrices were utilized as measures of uncertainty for the parameters.The norm accounts for covariances, while the trace takes into account only variances.Columns 5 and 6 of Table 1 show the norm and trace of the covariance matrix for each input function.As expected, one can observe that the uncertainty is generally lower for optimal input functions than for sub-optimal input functions.However, the relationship between the (a priori computed) D-optimality criterion and the (a posteriori determined) uncertainty is not monotonic; e.g., the uncertainty for the r = 5 optimal input function is slightly lower than that for the r = 6 input function while the D-optimality criterion value is greater for the r = 6 input function (Table 1).Furthermore, the uncertainty for the r = 6 sub-optimal input function shown is larger than that for the optimal r = 4 and r = 5 input functions even though the D-optimality criterion value is greater for the r = 6 sub-optimal input function.There are a number of possible reasons for this lack of monotonicity, such as nonlinearity of the IL-6 model and the fact that optimal experimental design theory is an approximate theory in the case of nonlinear models [8,10].
In order to corroborate that these findings also hold up when global rather than local analysis is used, the Morris method was implemented for computing sensitivity coefficients.This resulted in the optimal input functions listed in Table 2. Figure 2 shows the optimal input functions computed by both the local method and the Morris method for r from 1 to 6.It was observed that both methods resulted in optimal input functions with alternating IL-6 concentration levels, reminiscent of pseudo-random binary signals (PRBS signals).While the different methods sometimes result in input profiles that start high vs. low, the general shape of the functions is an oscillation between a low and a high value of the input.This suggests that PRBS-like signals may be favorable for minimizing parameter uncertainty in the IL-6 signaling model.An important take-away message from the simulation results is that optimal experimental design can be very effective at providing an a priori design, yet this does not replace evaluating the quality of a design a posteriori.

Discussion
The IL-6 model aims to capture the behavior of molecules involved in both the Jak-STAT pathway and the MAPK pathway [2].This model is larger and more nonlinear than most models that have previously been the subject of optimal experimental design involving signaling pathways [18,24,25].The present study shows that optimal experimental design can be utilized for improving the parameter accuracy of moderately complex nonlinear models.Further, the study shows that parameter uncertainty can be substantially decreased by optimizing the input function alone without simultaneously optimizing the experimental measurement time points for this particular system.This observation potentially has fortunate implications for experimental systems in which the input function can be controlled but the experimental measurement time points are dictated by practical considerations and cannot be changed.
The main step in the optimal experimental design was the sensitivity analysis for constructing the FIM.This involves constructing the Jacobian of the model equations so that the sensitivity equations can be integrated to obtain the sensitivity coefficients for the FIM.In principle, this method can be applied to any ODE model.The main step in the analysis is therefore to calculate all of the partial derivatives in the Jacobian, since the Jacobian will be different for every ODE model.For the Morris method, the information matrix is constructed by sampling from the parameter distribution and calculating finite difference derivative approximations for the sensitivities.In both cases-i.e., local and global methods-an ODE solver capable of accurately integrating the ODEs must be used, along with the optimization solver in order to evaluate the objective function for each optimization iteration.
Optimal experimental design using the D-optimality criterion was effective, according to the simulation results, in decreasing the parameter uncertainty in the model.However, one has to caution that this is only one aspect.For example, it has been shown in work by White et al. [36] that increased parameter accuracy does not guarantee an improvement in the accuracy of the model.
The form of the optimal input functions for the IL-6 model raises an interesting point regarding general system identification theory.Each of the determined optimal input functions takes the form of a PRBS-like sequence (see Figure 2), a commonly used input signal for system identification [37], even though this input function shape was not postulated during the formulation of the optimal design problem.Furthermore, computing an experimental design where sensitivities were computed via the Morris method also led to PRBS-like optimal inputs.Without commenting on generality, a PRBS signal seems to be a good choice for inputs of the investigated IL-6 signaling model, a result that may potentially carry over to other signaling pathway models.

Conclusions
Optimal experimental design was applied to the problem of minimizing parameter uncertainty in an IL-6 signaling model representing the Jak-STAT and MAPK pathways.The D-optimality criterion, operating on the FIM, was constructed using sensitivity equations, which were solved simultaneously with the equations of the model.Piecewise constant input functions were determined by solving this optimization problem; the piecewise nature of the inputs lays the groundwork for implementing the determined IL-6 concentration profiles on an experimental system.The optimal input functions resulted in decreased parameter uncertainty for the model, as observed from simulations in which parameters were fitted using the optimal input functions for inducing the signaling system.Interestingly, the determined optimal input functions took on the shape of PRBS signals even though this was not in any way postulated by the optimal experimental design problem.This observation was further validated by formulating and solving the problem using a global method in addition to the local method.Future work should corroborate these findings by applying the determined optimal experimental design to an actual experiment.
Table A1.Parameter values from initial parameter fit [2] and from optimal r = 5 input function.

Figure 1 .
Figure 1.Example of a piecewise constant input function for IL-6 (Interleukin-6) concentration with r = 4.The input function shown here can be represented as a vector of the four concentration levels in chronological order: c = [1.1,5.2, 0.9, 7.4] nM IL-6.

Figure 1 .
Figure 1.Example of a piecewise constant input function for IL-6 (Interleukin-6) concentration with r = 4.The input function shown here can be represented as a vector of the four concentration levels in chronological order: c = [1.1,5.2, 0.9, 7.4] nM IL-6.

Figure 2 .
Figure 2. Optimal input functions for values of r from 1 to 6. Blue indicates input functions identified as optimal from the local sensitivity method, green indicates input functions identified as optimal from the Morris method.(a) r = 1 (b) r = 2 (c) r = 3 (d) r = 4 (e) r = 5 (f) r = 6.

Figure 2 .
Figure 2. Optimal input functions for values of r from 1 to 6. Blue indicates input functions identified as optimal from the local sensitivity method, green indicates input functions identified as optimal from the Morris method.(a) r = 1 (b) r = 2 (c) r = 3 (d) r = 4 (e) r = 5 (f) r = 6.

Table 1 .
Optimal and sub-optimal input functions for values of r from 1 to 6.

Table 2 .
Optimal input functions for values of r from 1 to 6 via the Morris method

Table 2 .
Optimal input functions for values of r from 1 to 6 via the Morris method