Parameter Estimation in Computational Systems Biology Models: A Comparative Study of Initialization Methods in Global Optimization

— This paper compares different initialization methods and investigates their performance and effects on estimating kinetic parameters’ value in models of biological systems. Estimating parameters values is difficult and time-consuming process due to their highly nonlinear and huge number of kinetic parameters involved. Global optimization method based on an enhanced scatter search (ESS) algorithm is a suitable choice to address this issue. However, despite its resounding success, the performance of ESS may decrease in solving high dimension problem. In this work, several choices of initialization methods are compared and experimental results indicated that the algorithm is sensitive to the initial value of kinetic parameters. Statistical results revealed that uniformly distributed random number generator (RNG) and controlled randomization (CR) that being used in ESS may lead to poor algorithm performance. In addition, the different initialization methods also influenced model accuracy. Our proposed methodology shows that initialization based on opposition-based learning scheme have shown 10% better accuracy in term of cost function.


I. INTRODUCTION
Kinetic models of living cells have drawn the attention of both practitioners and researchers in recent years [1]. Their applications are important in metabolic and bioprocess engineering as they facilitate scholars to better understand, accurately predict and consistently improve the desired products in systems biology [2,3]. The models are formulated by means of ordinary differential equations (ODEs) to mimic various functional behaviours such as glycolysis reactions via metabolic pathway and phosphorylationin signal transduction of human cells. Due to the highly nonlinear biological systems, building such model is considered both challenges and time-consuming [4].
One important aspect of model building is parameter estimation, which consists of finding the best possible value of kinetic parameters that produce best fit model to the experimental data. The goodness of fit can be measured by minimizing distance value in the simulated model and experimental data. Thus, searching best parameter values in kinetic model can be depicted as a nonlinear optimization problem [5] and this class of problem is difficult to be solved. In this view, various optimization algorithms have been proposed in parameter estimation and their findings revealed that local optimization often fails to obtain snear-optimal solution [6]. Although improvements such as iterated local search have been proposed, they still consume high computational cost. Consequently, global optimization which is based on metaheuristic is an ideal option to address this issue. Global methods are quite capable in parameter estimation problem as they are more likely to reach the global minimum compared to local methods.
Enhanced scatter search (ESS) is one of the metaheuristic algorithm which have recently shown to yield promising outcomes in biological problems [7,8]. The algorithm benefits from global exploration and local exploitation using various choices of local search. The balanced tradeoff between global and local methods in ESS has shown promising results in solving optimization problems. However, when dealing with high dimension problem involving hundreds of kinetic parameters, performance of most global methods including ESS are deteriorate. One of the most neglected mechanisms in global methods is the way they generate the initial solution which were commonly derived using random number generator (RNG). The initialization methods may influence the efficiency and performance of the optimization algorithm in terms of its probability in finding the global minima, convergence's rate and variance of statistical results [9]. To date, only a few works have been done for comparing initialization methods in optimization. So far, no comparative study with regard to initialization method has been done in large-scale parameter estimation problem, particularly in the biological domain. The limitation of existing work in the field of global optimization is they only rely on RNG for initialization and only focus on search operator or the way new solution are produced. The high complexity of the problem such as in biological domain or healthcare is challenging and applying optimization method must properly select the best initialization because it will influence the www.ijacsa.thesai.org output. Hence, this issue motivates this research to further investigate the effect of different initialization method. This paper compares and investigates the effects of several initialization methods (also known as diversification generation method in ESS algorithm) from the context of parameter estimation in systems biology models. The evolutionary algorithm based on ESS is utilized in this study due to its efficiency and reliability in parameter estimation problem [7]. The paper is organized as follows: Section II explicates the problem statement in parameter estimations; Section III delineates ESS algorithm; Section IV introduces several initialization methods; Section V compares the methods and presents the discussion of their results and Section VI presents the conclusion of this study.

II. PROBLEM BACKGROUND
In a nonlinear kinetic model of biological systems, the parameter estimation problem deals with finding an unknown value of kinetic parameters to minimize a distance (objective or cost function) between simulated model and real data. The value of cost function determines the goodness of fit of the model. The observables, which are referred to as the output state variable, are experimentally measured. The cost function of this problem, which is also known as weighted nonlinear least squares is defined as: where is the number of experiments, is the number of observables per experiment and is the number of samples per observable in each experiment. Time series experimental data is denoted as and predicted model is denoted as . The kinetic parameters vector to be estimated is . The time span for observables is denoted as and finally represents the weight matrix to balance the contributions of the observables. Minimization of the above cost function is subject to the following constraints: where derivative of ̇ is the function of ODEs model that describes the dynamics of biological systems, is the initial condition at time , is the observable functions and and are the lower bound and upper bound of the kinetic parameter vector respectively. This nonlinear and multimodal problem consists of many local minima. Thus, the process of finding the global minima is both challenging and time-consuming.

III. ENHANCED SCATTER SEARCH (ESS) ALGORITHM
An enhanced scatter search (ESS) is a metaheuristic that belongs to the family of evolutionary algorithms. This algorithm is similar with genetic algorithm (GA) with regards to maintaining and updating their population members and evaluating their cost function in an iterative cycle. However, unlike GA, ESS does not use crossover and mutation as their evolutionary operators. Instead, it uses the combination among members in a reference set (RefSet). In this study, four phases of ESS algorithm are used, namely: 1) initialization method, 2) RefSet update method, 3) RefSet member generation and combination method, and 4) hybrid of the local search method. More advance designs and their mechanism can be found in [10,11].
This algorithm starts with randomly generating population of diverse vectors by means of initialization (diversification generation) method. The size is ten times the problem size to ensure that the large initial solutions in the search space are widely sampled, thus increasing the chances of avoiding local minima. Although uniformly distributed random number is a popular method usually utilized to generate initial solutions in various optimization algorithms, there are other strategies that may provide better initial solutions. Therefore, we compared and investigate four different initialization methods which will be briefly discussed in Section IV.
After the diverse vectors are generated, each vector is evaluated and half of the RefSet members ⁄ ( is the RefSet size) is formed. The diversification method produces high quality initial RefSet member. The remaining RefSet members are choosen from the RefSet by random cycle to complete a RefSet. Then, the subset generation produces pairs of members in RefSet. Let us consider members of a RefSet, , to be combined with the rest of members in RefSet, [ ] . The pairs of the combination ( and ) are defined as follows: and | | Every pair of combination ( and ) in the RefSet members is used to create new hyper-rectangles which are defined by their relative positions and distance and thus, resulting in a new solution within them. The hyper-rectangles based combination methods are applied and are defined in the following equation: (11) where is new solution generated and is the random number, [ ] . This combination strategy is similar to the mutation operator in differential evolution (DE) [12,13], which is effective in updating population members. In ESS, www.ijacsa.thesai.org the vectors of combination ( and ) are systematically generated. They are not randomly generated, as practiced in DE. Using this combination strategy, every RefSet member generates a hyper-rectangle among the rest of RefSet member. The new number of solution produces solution for each RefSet member. After this strategy is implemented, a new solution (offspring) is generated with different distance and direction around their RefSet members (parents). In this case, if the offsprings have better (lesser) fitness value compared to their parents, the current solutions will be replaced. Otherwise, the same RefSet members will be used for the next iteration. In order to accelerate convergence, gradient local search is performed using Sequential Quadratic Programming (SQP). The algorithm is applied using fmincon solver in MATLAB. This solver minimizes the cost function using the results obtained in ESS using different vectors. If the solution obtained by fmincon outperforms the solution generated by ESS, the solution from fmincon will replace the current solution and it will in turn be added to RefSet members for further update. Otherwise, the solution from fmincon will be discarded. This process is repeated until the stopping criteria are met.

IV. INITIALIZATION METHODS
We implement five initialization methods in the ESS algorithm in order to compare and investigate their effects on parameter estimation. The methods are random number generator (RNG), controlled randomization (CR), oppositionbased learning (OBL), quasi-opposition learning (QOBL) and chaotic (Tent) map.

A. Random Number Generator (RNG)
The most commonly used initialization method in optimization algorithms is random number generator (RNG). RNG is defined as below: Let ) be the th member of the population, each is generated between lower and upper bound ( , ). In summary, it generates uniformly distributed random numbers as in the following equation: where is the random numbers between 0 and 1. The vector of contains a list of random initial population generated between lower bound and upper bound [ , ] for each of the variable.

B. Controlled Randomization (CR)
Unlike RNG, controlled randomization (CR) strategy generates the first five populations ( ) of equal size for each vector as in the following equation [14]: (13) where is the vector of candidate solutions, is the random numbers in the between 0 and1, and is the number of kinetic parameters. After the first five vectors are generated, the remaining vectors are generated randomly and all initial solutions are put in the boundaries: (14) where and are lower and upper bounds, respectively. This strategy generates a set of diverse vectors which contain equal sizes of range in the first five vectors and other random vectors lie in sixth vector to diverse vectors. It should be noted that ESS algorithm used CR strategy as its default initialization method [11].

C. Opposition-based Learning (OBL)
Opposition-based learning (OBL) is introduced in the field of computational intelligence [15]. This scheme is subsequently applied in optimization areas [16]. The basic idea of OBL is to generate a set of opposite numbers from first initial solutions generated by RNG, as follows: Let  [ ] is a random value. The opposition value of is defined by: Both and ̆ is merged. Now, let us assume is the cost function in minimization problem. If cost function value of ̆ is smaller than , ̆ , point can be replaced with ̆. Otherwise, point will stay in the current population. All the population members will be evaluated and the initial population with fittest members among and ̆ is formed.

D. Quasi Opposition-based Learning (QOBL)
Another family of OBL is quasi-opposition learning (QOBL) which modified version of OBL that increases population uniformity [17]. Considering the opposite point in equation 16, the middle point { } is calculated as follows: Then, quasi-opposite point ̆ is selected randomly within the range of opposite points of ̆ and middle point : where is a random value drawn uniformly in the range of lower bound and upper bound. Like OBL, ̆ and are merged and the best solution is chosen, that is, the fitness of the cost function.

E. Chaotic Map
Another alternative of uniformly distributed random numbers for diversification generation is chaotic map. This approach is based on the deterministic and chaotic systems and it is not necessarily random. In this paper, we investigate (IJACSA) International Journal of Advanced Computer Science and Applications, Vol. 13, No. 8, 2022 476 | P a g e www.ijacsa.thesai.org one family of the chaotic map, which is Tent map [18] which is defined as: where is th variable of th individual in th iteration. is the initial variable that is generated randomly using RNG. In this strategy, solutions which are generated from Tent map are not predictable and are highly sensitive to initial variables.

V. RESULT AND DISCUSSION
A large-scale model is used to test the different initialization methods in ESS algorithm. The model involves dynamic processes that reproduce the response to a pulse in extracellular glucose concentrations of central carbon metabolism (CCM) in E. coli. This model consists of 18 metabolites: 17 internal metabolites in cytosol and 1 extracellular metabolite (glucose) in extracellular compartment.
These metabolites consists of PEP, G6P, PYR, F6P, G1P, 6PG, FDP, GAP, CPEP, CG6P, CPYR, CF6P,GLCex, CG1P, CPG, CFDP, CGAP and Glucose. The model also contains 48 reactions coupled with 166 kinetic parameters. The mathematical formulation and description of this model can be found in [19]. Table I   In order to obtain statistically significant result, we ran each initialization method discussed in Section IV, 20 times and reported the best, mean, and worst results; as well as average function evaluations, CPU time and standard deviation. Function evaluation for each run was limited to 100,000 (the stopping criteria) to let the algorithm obtain the best parameter values. The RefSet size used was 36, which is the recommended size in this problem. With the high number of function evaluations and hundreds of parameters, the minimization process is expected to consume very lengthy CPU time. To surmount this drawback, Parallel Computing Toolbox in MATLAB has been used and it expedited the computation by assigning each run to eight different processors (logical cores) simultaneously. In this strategy, eight computations for each method was run independently using parfor loop which is available from the abovementioned toolbox. It should be noted that a single run takes approximately 11 hours, so 20 runs take approximately 220 hour. Using the parallel strategy, 20 runs only take approximately 33 hour, which reduced 72.6% of CPU time needed. All methods were experimented on i7 CPU with 16GB RAM which implemented in MATLAB 2015. Table II shows that the best (minimum) cost function was obtained from QOBL method with . The second best value is 229.1855, which was obtained from CR. Only RNG, OBL and TENT produced cost function values which were slightly higher than the published benchmark value, 233.90. The results revealed that QOBL is the best method in finding global minimum. However, although QOBL presented the minimum value, its average standard deviation was relatively higher than RNG, OBL and TENT. RNG is the most consistent method followed by OBL, having 4.5955 and 4.7331 standard deviation each, respectively. In terms of search effort, RNG produced the lowest average of function evaluations (1.2327e+05) and also its CPU time is also the lowest with 3.8191e+04 seconds. It should be observed that QOBL is the best initialization method if we consider its ability in minimizing the cost function in largescale parameter estimation problem. Additional information to compare the different initialization methods is given in Fig. 1 and Fig. 2. The figures depicts the best curves (with minimum cost function ) among the 20 runs obtained from RNG, CR, OBL, QOBL and Tent methods. The curves show that all methods are able to minimize the cost function at a similar rate in terms of function evaluations and CPU time. Note that we have set the same initial guess as the initial value for all methods. This gives fairer comparison and assumes that the search space is feasible. In Fig. 1, starting from the first fractions of 1,000 function evaluations, QOBL has better speed and found acceptable cost function value when it reached approximately 25,000 function evaluations. Meanwhile, CR has the slowest speed until it reached around 40,000 evaluations. All methods continued to progress they reached the final evaluations. In this case, QOBL found the best value of 210.051 at 123,979 function evaluations. At the end of the evaluations, all solutions were able to achieve equivalent solutions in terms of quality, although, the best (minimum) value was obtained by QOBL while the worst (maximum) value obtained by RNG. In www.ijacsa.thesai.org Fig. 2, default initialization method based on CR has obtained slow convergene rate compared to others. It can be noticed where CR obtained acceptable value of cost function when CPU time reached nearly 8 hours, while QOBL reached the acceptable value at nearly 1.3 hours. The results revealed that RNG which is mostly used initialization method and CR as default initialization method in ESS have obtained poor result. Meanwhile, initialization method based on QOBL is a better alternative which was not only able to speed up convergence, but also obtain optimal solution.
To evaluate the quality of the parameter estimates, we compared the best result (QOBL) with a published benchmark result [20]. One thing to note is benchmark result used conventional CR as their initialization methods for the ESS. Table III shows that our study produced the best cost function, = 210.05 compared to the benchmark's = 233.90. However, our work obtained a bigger number of function evaluations compared to CR, with a difference of 33,251 evaluations. In terms of CPU time, the benchmark also produced shorter time of approximately 3 hours for a single run. Due to different stopping criteria used (the benchmark study uses CPU time) and different hardware specifications, comparing QOBL and CR in terms of CPU time seems unfair. It should be noted that in terms of efficiency in finding global minimum, method used in this study produced better results compared to CR.   The decision variables (kinetic parameters) obtained from QOBL may provide the optimal model prediction since it produced the lowest cost function and may produce best fit to experimental data. The goodness of fit can be measured by calculating root-mean-square-error ( ) for all metabolites.
is used to measure prediction error which is the different between experimental data and predicted model. The following equation defines , with the same notation defined in Eq. (1). In this case, normalized RMSE is used to cater for different magnitudes of observables. Each RMSE is divided by the range of value of observables as defined as,  Table II shows our is lower than the benchmark, means that it has better fit compared to the parameters obtained in the benchmark.
To measure the goodness of fit from parameter estimates with QOBL, we plot the model prediction over experimental data as shown in Fig. 3. For the sake of brevity, we plot only four out of nine metabolites. We chose metabolites which have very high nonlinear biological system, namely, PEP, G6P, PYR and F6P. The figure shows dynamic concentration change of extracellular glucose that responded to a pulse in central carbon metabolism. To compare the goodness of fit with other parameters, we also plot another fit based on parameters value from original published results and benchmark results. It should be observed that kinetic parameters retrieved from initialization methods based on QOBL represent a good fit between experimental data and predicted model.

VI. CONCLUSION
This paper studies different initialization methods and investigated their performance and effects on solving largescale parameter estimation problem. We compared five initialization methods which are based on stochastic and randomization methods and implement them in ESS algorithm. Experimental results revealed that the choice of initialization (diversification generation) methods influenced the performance of the algorithm. The quality of solution, speed of convergence and statistical results were obtained with different characteristics derived from different initialization methods. Our statistical analyses revealed that the most popular initialization method, random number generator (RNG) performs poorly and there are significant better alternatives to this method, which have comparable computational requirements. In addition, the accuracy of model prediction also depends on the choice of initialization methods. Further investigation is needed to discover whether the same findings can be produced when different models and problems associated with parameter estimation are used. More intensive studies also need to be conducted on why some methods are generated more consistent performance in terms of statistical analysis and value of kinetic parameters in biological systems.