**1. Introduction**

In silico optimization is a rapid and cost-effective method for finding optimal solutions in optimization problems. The two currently available methods are local and global methods. Each, however, comes with inherent systematic problems that require troubleshooting. The local optimization method performs poorly when solving non-linear and dynamic biological problems [1]. In this type of problem, the global search is a more effective method as it is capable of finding an optimal solution that satisfies all requirements. In fact, the global search method has recently been subject of much attention within the scientific community [1–3]. Parameter estimation is an essential phase in the simulation of the biological system because the value of the parameter determines the behaviour of the model. In estimating parameter values, it is important to first determine the objective function that can then be minimized by using suitable optimization methods. The objective function must be minimized to obtain the ideal value of the parameter. Several optimization methods have, therefore, been applied to estimate parameter values [4]. After comparing several methods, Baker et al. [4] pointed out that while the Genetic Algorithm (GA) was less time-consuming than Simulated Annealing (SA), it faced the local minima problem. In MCMC approaches are often applied to determine the posterior distribution of rate parameters. However, developing a good MCMC sampler for its multimodal and dimensional parameter distribution is challenging. Valderrama-Bahamóndez and Fröhlich [5] found that parallel adaptive MCMC performed better in parameter estimations after comparing the ability of di fferent MCMC approaches to estimate the kinetic rate parameters of ordinary di fferential equation (ODE) systems.

In addition, global optimization algorithms, like Particle Swarm Optimization (PSO) [6,7], SA [8], and Scatter Search [9,10], have been successful in the estimation of parameters in di fferent biological models. In a comparative study on computational intelligence methods, carried out by Tangherloni et al. [11], the performance of various meta-heuristics was compared; in particular, their ability to estimate the parameters with a set of benchmark functions and synthetic biochemical models. Tangherloni et al. [11] concluded that classic benchmark functions lacked the comprehensibility that compounded the real-world optimization problem. The hybridization of optimization methods has been proposed so that the best features of each method could be utilized. Convergence to global optima is a consistent problem in standard PSO for multimodal and high-dimensional functions [12]. Hence, Fu et al. [12] proposed a hybrid method by incorporating the evolutionary operations of the Di fferential Evolution (DE) algorithm to improve its conventional velocity updating strategy in PSO. Harmony Search (HS) is a simple method with an e fficient and evolutionary algorithm. It has parameters like the Harmony Memory Considering Rate (HMCR) and Pitch Adjusting Rate (PAR) that solve the local optima problem. Furthermore, HS is able to balance between intensification and diversification whereby promising regions and non-explored regions are thoroughly and evenly explored in both processes [13]. This paper introduces a hybrid of Particle Swarm Optimization [14] and Harmony Search [15] (PSOHS) and simulates the essential amino acid metabolism for kinetic parameter estimation.

PSO estimates kinetic parameter values by using swarm intelligence. The candidate solutions are analogous to particles flying with specific velocities in specific directions in the search space, either alone or together with their companions. Presumably, the particle will move to its best positions. HS, on the other hand, is based on the analogy with natural musical performance processes. By improving the pitches of the instrument of each music player, a better harmony is achieved and, hence, the quality of the performance [16]. Gao et al. [16] further discussed this method and its applications. The historical development of the HS algorithm structure had been extensively reviewed by Zhang and Geem [17].

This paper introduces a new hybrid algorithm, PSOHS, into the SBToolbox to estimate kinetic parameters by simulating the aspartate metabolism of isoleucine, lysine, and threonine in a small plant of the mustard family, *Arabidopsis thaliana*. Its model is often chosen in research because it has been demonstrated to render better results when analysing plant development and growth. This paper focuses on the biochemical reactions of essential amino acid production in *Arabidopsis thaliana*.

#### **2. Materials and Methods**

This section describes and discusses the problem formulation, and it details the hybridization of PSO and HS. A basic PSO is also described to di fferentiate between the basic PSO and the hybrid version.

#### *2.1. Problem Formulation*

This paper uses a biochemical system, which is a mathematical modelling framework based on ordinary di fferential equation (ODE). In the system, biochemical processes are represented using power-law expansions in the variables. In a biochemical process, a system of kinetic equations is

formed according to the information regarding the underlying network structure of a pathway, and the parameters are acquired from the literature or estimated values are acquired from a data fit. Chemical kinetics mean rates of chemical reactions. The parameter estimation problem is posed in Equation (1) below. There, *s(X)* denotes a biological compound, depending on a set of parameters *X* = *(X1,X2,X3,* ... *Xd)*, where *d* is the total number of parameters. Hence, the reaction rate of the compound s is presented by

$$\begin{array}{l} \frac{ds}{dt} = \mathcal{g}(s(X), t), \\ s(t\_0) = s(0), \\ y = \mathcal{g}(s(X), t) + \varepsilon. \end{array} \tag{1}$$

In Equation (1), *g* represents a nonlinear function and *t* represents the sampling time. Then, *y* is a time series of simulated data, also known as the output of the model, and *e* represents the noise data that is randomly produced by Gaussian noise *n* (1,0). The purpose of the parameter estimation is to discover the set of the optimal parameter, which is denoted as *X*. Then the variance between the simulated time-series data denoted as *y* and experimental time-series data denoted as *yexp* can be reduced. The variance is calculated by applying the nonlinear least squared error function, *f(X)*, which is shown in Equation (2):

$$f(X) = \min \sum\_{i=1}^{n} (y^{exp} - y)^2 \tag{2}$$

In Equation (2), *n* is the total number that maximum value generated and *i* is the index variable.

#### *2.2. Particle Swarm Optimization (PSO)*

PSO is a swarm intelligence-based optimization algorithm. The method was developed by Kennedy and Eberhart [6] and had been inspired by the social behaviour of flocking birds and schooling fish when searching for food [6].

In PSO, particles act as a possible solution in the search space of a problem. Every particle in the search space is assigned a certain velocity so that its movements through the search space can be determined. In the search space, each particle's movement is a ffected not only by its local best-known position, but it is also directed toward the best-known positions, where the best conditions are those that have been found by other particles.

Hence, there are two positions in the whole swarm, the local best-known position and the best-known position in the swarm. *pbest* represents the local best-known position while *gbest* represents the best-known positions in the swarm.

Figure 1 shows the PSO flowchart. First, each particle is assigned a random position within the problem space to form an initial population. Then, the fitness of each particle is evaluated and compared with *pbest*. If the current value is better than *pbest*, then this value is updated to the new *pbest*. Then, the particle's best-known position and the swarm's best-known position are updated. The termination criterion is when the maximum number of iterations is performed, or a solution meets the adequate objective function value.

**Figure 1.** Flowchart of Particle Swarm Optimization (PSO).

#### *2.3. Harmony Search (HS)*

Geem et al. [13] developed HS after being inspired by the music improvisation process. The algorithm achieved better harmony the same way as music players improved the pitches of their instruments [18]. To escape from the local minima, HS performs a stochastic random search instead

of a gradient search. There are three rules to generate perfect harmony, which are Random Selection (RS), Harmony Memory Considering Rate (HMCR), and Pitching Adjust Rate (PAR). Figure 2 shows the flowchart of HS. First, a population of random harmonies in Harmony Memory (HM) is initialized. Several randomly generated solutions are included in the initial HM. The next step is to improvise a new solution from the HM followed by updating the latter. In each repetition, the algorithm improvises a new harmony, and the latest harmony is generated by the following three rules: (1) HMCR is used to select the variables of the new harmonies from the overall HM harmonies; (2) PA is responsible for local improvement; and (3) RS provides random elements for the new harmony. If the latest harmony is better than the existing one, then it evaluates and replaces the worst harmony in HM. This process is iterated until the stopping criteria are met.

**Figure 2.** Flowchart of Harmony Search (HS).

#### *2.4. A Hybrid of PSO and HS (PSOHS)*

PSOHS is proposed in this paper, which is a hybrid of HS and PSO. Figure 3 shows the steps taken by the proposed algorithm in order to acquire the optimal parameter values.

**Figure 3.** Steps in PSOHS for estimating parameter value. The red box is Harmony Search, which is hybridized with PSO to improve its performance.
