what value of s as a fraction of km is required to obtain 20 vmax

Comput Methods Programs Biomed. Author manuscript; available in PMC 2015 Feb 1.

Published in last edited form every bit:

PMCID: PMC3906649

NIHMSID: NIHMS533785

Statistical Identifiability and Convergence Evaluation for Nonlinear Pharmacokinetic Models with Particle Swarm Optimization

Seongho Kim

¹Biostatistics Core, Karmanos Cancer Plant, Wayne State Academy, Detroit, MI, 48201

Lang Li

^twoDepartment of Medical and Molecular Genetics, Indiana University School of Medicine, Indianapolis, IN, 46032

Abstract

The statistical identifiability of nonlinear pharmacokinetic (PK) models with the Michaelis-Menten (MM) kinetic equation is considered using a global optimization approach, which is particle swarm optimization (PSO). If a model is statistically not-identifiable, the conventional derivative-based estimation arroyo is often terminated earlier without converging, due to the singularity. To circumvent this difficulty, we develop a derivative-complimentary global optimization algorithm by combining PSO with a derivative-gratis local optimization algorithm to meliorate the rate of convergence of PSO. We further propose an efficient approach to non merely checking the convergence of estimation only besides detecting the identifiability of nonlinear PK models. PK simulation studies demonstrate that the convergence and identifiablity of the PK model can be detected efficiently through the proposed approach. The proposed approach is then practical to clinical PK data forth with a two-compartmental model.

Keywords: Michaelis-Menten kinetic equation, nonlinear models, particle swarm optimization, pharmacokinetics, statistical identifiability

1. Introduction

The nonlinear modeling is a routine simply absolutely necessary statistical method in analyzing drug concentration information measured over time in pharmacokinetics (PK). In PK studies, Michaelis-Menten (MM) equation is often employed to draw the intrinsic clearance

where Vmax is the maximum enzyme activity; Km is an changed function of the affinity betwixt drug and enzyme; C(t) is an unbound drug concentration. Km is also called the MM abiding having the units of C(t). The deterministic and statistical identifiabilities of parameters in the MM equation have been examined (Tong and Metzler, 1980; Metzler and Tong, 1981; Godfrey and Fitch, 1984). The deterministic identifiablity is concerned with whether the model parameters can be identified with racket-costless data, while the statistical identifiability is the possibility of identifying the model parameters with noise data.

Although numerous methods have been presented to detect the non-identifiable parameters deterministically, such equally the Laplace transform (Godfrey and DiStefano, 1987), the similarity transformation approach (Vajda et al., 1989), the Voterra and generating power serial approaches (Lecourtier et al., 1987), the differential algebra approach (Saccomani et al., 2003), and the alternate conditional expectation algorithm (Hengl et al., 2007), there has been much less evolution in statistical identifiability analysis of PK models. 1 of the empirical approaches to assessing the statistical identifiabilty is the local sensitivity assay. The local sensitivity analysis in the statistical identification uses the first partial derivatives of the differential equations with respect to the parameters, and depends on the not-singularity of the Fisher information matrix, which is equivalent to the Taylor serial method and differential algebra method (Hidalgo and Ayesa, 2001; Wynn and Parkin, 2001).

Still, the local sensitivity analysis is probable to make a wrong decision if the estimate is far from the truthful value or the model has very complicated dynamics. Yue et al. (2008) thus proposed the global sensitivity analysis for robust experimental pattern based on the modified Morris method (Morris, 1991), but it notwithstanding requires an initial guess or prior knowledge concerning the underlying relation of the parameters. Therefore, we advise an arroyo not only to accessing the identifiability globally only also to requiring no preprocessing to obtain an initial estimate or prior knowledge.

A number of estimation approaches were developed for population PK analysis (Beal and Sheiner 1982; Lindstrom and Bates 1990; Vonesh and Carter 1992; Wolfinger, 1993; Kim and Li, 2011). Nearly approaches are a derivative-based local optimization method, however. A well-known challenge of the local optimization, such as the Newton and akin methods, is stuck at the saddle points or a local optimum so that the initial values are required to lie inside a relatively small neighborhood of the true optimum to notice a global optimum, and the derivative-based method is often terminated earlier due to the singularity. The singularity trouble tin become more prominent when the model is statistically not-identifiable. These issues urge us to use a derivative-free global optimization algorithm since information technology tin can avert the singularity problem besides as seek the best parameter estimates of nonlinear models regardless of the presence of multiple local optima.

One interesting evolution based global optimization approach, particle swam optimization (PSO), was developed by Kennedy and Eberhart (Eberhart and Kennedy, 1995; Kennedy and Eberhart, 1995). PSO algorithm is a derivative-free approach and becoming very popular due to its simplicity of implementation and robust convergence capability. Using PSO algorithm, Kim and Li (2011) developed a global search algorithm, P-NONMEM, for nonlinear mixed-effects models to meet the challenges of the local optimization in NONMEM, which is one of the near popular approaches in PK studies. However, NONMEM uses a Broyden-Fletcher-Goldfarb-Shanno (BFGS) quasi-Newton algorithm, which is a derivative-based approach, so that it is non complimentary from the singularity problem. For this reason, nosotros develop a modified version of PSO algorithm, which is the PSO coupled with a derivative-free local optimization algorithm (LPSO), in lodge to estimate the parameters regardless of the identifibility.

1 claiming of PSO algorithm is the lack of convergence criteria. The number of function evaluations is often used as a stopping benchmark forth with incorporating the choice of a problem-dependent parameter, which relies on the gradient or divergence between the previous and the current estimates. However, this arroyo doesn't accept the random or stochastic behavior of PSO into account so that it will make the estimation stopped before reaching a global optimum. It as well focuses but on the identifiable situations. Therefore, it is desirable to have a reliable convergence criterion for detecting when the optimization process has found the global optimum even for not-identifiable conditions. Nosotros thus propose several approaches to non only diagnosing the convergence of PSO but also detecting the statistical identifiability.

In Section 2, a brief description of a two-compartment model with Michaelis-Menten kinetic equation is given. The nonlinear PK models with PSO are introduced in Section 3. In Department iv, the proposed PSO algorithm and its convergence criteria are described in details. Simulation studies are performed to evaluate the proposed approaches and real clinical PK data and so are practical in Section 5. In Department 6, conclusions are reached.

2. Michaelis-Menten Kinetic Equation and Ii Compartmental Pharmacokinetics Model

Statistical Identifiability with the Michaelis-Menten kinetic equation

Information technology is well known that the drug metabolism rate follows the Michaelis-Menten (MM) kinetics equation:

$V (t) = \frac{d C (t)}{d t} = \frac{51000 a x \cdot C (t)}{Yard m + C (t)},$

where V(t) is the velocity of the reaction, Vmax is the maximum velocity, Km is the MM constant, and C(t) is the drug concentration. Monod (1949) first applied the MM equation to microbiology for the growth rate of microorganisms.

The MM equation by and large describes the relationship between the rates of substrate conversion by an enzyme to the concentration of the substrate. In this relationship, V(t) is the charge per unit of conversion, Vmax is the maximum rate of conversion, and C(t) is the substrate concentration. The MM abiding Km is equivalent to the substrate concentration at which the rate of conversion is half of Vmax. Km approximates the affinity of enzyme for the substrate. A small Km indicates high affinity, and a substrate with a smaller Km will approach Vmax more quickly. Very high C(t) values are required to approach Vmax, which is reached but when C(t) is high enough to saturate the enzyme (Hein and Niemann, 1962).

In pharmacology research, the statistical identifiablity often occurs with the MM equation. Suppose the observed data y(t) follows a normal distribution with the MM equation at a fourth dimension indicate t given the parameter θ =(Vmax, Km):

y(t) ∼N D{l o chiliad f(θ,t),σ ^two},

where f(θ, t) = V(t) and ND stands for a normal distribution. Withal, when Km is much higher than the concentration C(t) (i.e., Km ≫ C(t)), the office f(θ, t) is shut to $\frac{Five thousand a x}{K yard} \cdot C (t)$ in the equation below:

$f (θ, t) = \frac{5 m a x \cdot C (t)}{K grand + C (t)} \approx \frac{Five m a x}{Thou m} \cdot C (t) i f C (t) ≪ K m .$

In add-on, when Km is much smaller than the concentration C(t) (i.e., Km ≪ C(t)), f(θ, t) is shut to Vmax in the equation beneath:

$f (θ, t) = \frac{V m a x \cdot C (t)}{K chiliad + C (t)} \approx V m a x i f C (t) ≫ K m .$

In other words, if the concentration C(t) is much either less or greater than Km, i will non be able to estimate both Km and Vmax separately due to identifiability.

Ii Compartmental Intravenous Pharmacokinetic Models with the Michaelis-Menten kinetic equation

Compartmental PK assay uses kinetic models to describe and predict the concentration-time curve for both oral (PO) and intravenous (4) administration. PK compartmental models are often similar to kinetic models used in other scientific disciplines such as chemical kinetics and thermodynamics. The simplest PK compartmental model is the ane-compartmental PK model with oral dose administration and first-order elimination (Chang, 2010). A two-compartmental IV model with the MM equation is considered for this study. In this case, its PK is described by the system of the ordinary differential equations (ODEs):

$\begin{matrix} \frac{d A_{1} (t)}{d t} = - C L \cdot \frac{A_{1} (t)}{V_{1}} + C L_{12} \cdot (\frac{A_{2} (t)}{V_{2}} - \frac{A_{i} (t)}{5_{1}}); \\ \frac{d A_{two} (t)}{d t} = - C {Fifty}_{12} \cdot (\frac{A_{two} (t)}{V_{ii}} - \frac{A_{1} (t)}{V_{i}}); \\ C Fifty = \frac{Q h \cdot C Fifty i n t}{Q h + C L i n t}; \\ C L_{i n t} = \frac{V m a x}{One thousand k + \frac{A_{one} (t)}{V_{1}}}; \\ (A_{1} (t), A_{ii} (t)) ∣_{t = 0} = (D o s e, 0), \end{matrix}$

where (A ₁(t), A ₂(t)) are amounts of drug in systemic and peripheral compartments at time t, respectively, (V ₁, 5 _two) are volumes of distribution in systemic and peripheral compartments, respectively, CL ₁₂ is the inter-compartment rate constant, CL is the systemic clearance, CL _int is the intrinsic hepatic clearance, Vmax is the maximum of velocity, Km is MM constant, and Qh is the hepatic blood menstruation known every bit 80 50/h.

Because the ODEs are nonlinear, there exists no closed-form solution and a numerical approach should be used to solve the differential equations. We utilise the R package odesolve to deal with the ODEs. Due to the nature of the clinical written report, but the systemic concentrations are observable from PK study and its predicted concentration at time t is given past

where θ = (logV ₁, logV ₂, logCL ₁₂, logVmax, logKm).

iii. Pharmacokinetic Nonlinear Models and Particle Swarm Optimization

We illustrate the nonlinear PK model with the MM equation in this department. The observed drug concentration is described by a non-linear model,

l o 1000 y _i ∼N D(50 o g f(θ,t _i),σ ^two), i = 1,…,Northward,

where North is the number of time points, y_i the drug concentration at time t_i , and f(·) a nonlinear part of population PK parameter vector θ, and its log-transformed value $l o g f (θ, t) = l o g \frac{A_{1} (t)}{V_{1}}$ . The PK model is oft assumed to follow a log-normal mistake model since the observed systemic concentration y_i is greater than zero. And then the log-likelihood function for (θ, σ ²) is

$- \frac{N}{ii} l o m (2 π) - Due north l o 1000 (σ^{2}) - \frac{ane}{two σ^{2}} \sum_{i = i}^{N} {(fifty o one thousand y_{i} - l o g f (θ, t_{i}))}^{ii} .$

Since the function f(·) has no airtight form solution, the parameters are estimated using numerical methods, such as Newton-type approximation, Laplace asymptotic approximation, and Markov concatenation Monte Carlo simulation. Yet, these approaches are derivative-based methods so that it might not be complimentary from singularity, especially when the model is non identifiable. Furthermore, most derivative-based algorithms are local optimization approaches. For these reasons, we accommodate a global derivative-free optimization algorithm, particle swarm optimization (PSO), to deal with both identifiable and non-identifiable models.

PSO was originally developed by Kennedy and Eberhart (1995) as a population-based global optimization method. Its evolutionary algorithm stochastically evolves a group of particles. PSO allows each particle to maintain a retentiveness of its all-time fitting. Each particle's trace in the search space is and then determined by its ain retentivity of best fittings. Individual particle moves towards a stochastically weighted average of these positions, until they converge to the global best. It is used to solve a wide array of dissimilar optimization problems because of its attractive advantages, such as the ease of implementation and its gradient free stochastic algorithm. It has been proved to be an efficient method for many global optimization bug, and not suffering from the difficulties encountered by other evolutionary computation techniques. For instance, PSO does non endure from some of genetic algorithm (GA)'south difficulties, such equally interaction in the grouping enhances rather than detracts from progress toward the solution. In addition, PSO has retentiveness, which GA does not accept. Change in genetic populations results in destruction of previous noesis of the problem. In PSO, individuals who fly past optima are tugged to return toward them, pregnant that knowledge of good solutions is retained by all particles (Kennedy and Eberhart, 1995; Eberhart and Kennedy, 1995). For an overview of PSO and its variants, see Englbrecht (2007).

Allow particle s exist an chemical element of the population. Its position vector is $x^{s} = {(x_{g}^{s})}_{k = ane, \dots, Grand}^{'}$ , and updating velocity vector is ${five}^{s} = {(v_{k}^{due south})}_{yard = i, \dots, M}^{'}$ , where K is the total number of iterations of PSO and S is the population size (s = ane, ..., S). Its best previous positions of itself (i.e. local best) and the population (i.e. global best) are represented as $x_{50 b eastward southward t}^{due south}$ and x_gbest , respectively. The velocity $v_{k + 1}^{s}$ and the position $10_{chiliad + one}^{s}$ at the (m+1)thursday iteration of particle s are calculated according to the following equations:

$5_{1000 + i}^{s} = w_{g} 5_{k}^{due south} + c_{i} r_{1} (x_{50 b e south t}^{s} - 10_{k}^{southward}) + c_{21} r_{2} (x_{g b e due south t} - 10_{k}^{s});$

(4)

where due west_yard is called inertia weight (0 ≤ due west_thousand ≤ 1), c ₁ and c ₂ are the two positive constants chosen cerebral and social coefficient, respectively, r _one and r _ii are the two random sequences in the range [0,i], and k is the iteration number. The low values of constants c _i and c ₂ allow each particle to roam far from the target regions before beingness tugged back, while the high values effect in abrupt movement towards target regions. For this reason, these constants are conventionally set as 2.0, equally we did in our simulation studies. The inertia weight w_k is

$w_{m} = w_{chiliad a x} - \frac{k}{K} ({westward}_{one thousand a x} - w_{1000 i northward}),$

(6)

where due west_min and w_max are user-defined constants in the range [0,1] and w₁₀₀₀ = west_min ≤ due west_k ≤ west_max = w ₀. The inertia weight is adapted to control the impact of the previous history of velocities on the current velocity and to influence the trade-off between global (wide-ranging) and local (nearby) exploration abilities of the "flying points". A larger inertia weight facilitates global exploration (searching new areas) while a smaller inertia weight tends to facilitate local exploration to fine-tune the current search area. Suitable selection of the inertia weight tin can provide a balance between global and local exploration abilities and thus crave less iteration on average to find the global optimum (Shi and Eberhart, 1998). In order to use these properties of the inertia weight, we use a dynamic inertia weight by linearly decreasing it equally described in (6) and so that PSO tin can escape from premature convergence when it gets stagnated (Zhang and Cai, 2009).

4. LPSO: Particle Swarm Optimization Coupled with a Local Optimization Algorithm

We advise a modified PSO coupled with a local optimization algorithm to better the rate of convergence and call the proposed algorithm LPSO. Kim and Li (2011) proposed P-NONMEM that is a combined approach between PSO and NONMEM for mixed-effects models. Notwithstanding, NONMEM, one of the most popular algorithms for PK study, is a derivative-based algorithm then that information technology will not avoid the singularity trouble of non-identifiable models. For this reason, we comprise a derivative-free local optimization algorithm, Nelder-Mead method, into PSO to deal with the local all-time. Information technology is a derivative-free direct search method based on evaluating an interesting office at the vertices of a simplex iteratively past shrinking the simplex to find meliorate points until some desired leap (Nelder and Mead, 1965). It is also called a simplex search algorithm adult past Nelder and Mead (1965). Note that the term simplex is a generalized triangle in a certain dimension. Nelder-Mead method requires no derivative data, making it suitable for problems with not-smooth functions or/and discontinuous functions. Its full general algorithm is composed of the following two steps: construct the initial working simplex and repeat the transformation of the working simplex until it converges. At that place are four transformations to compute the new working simplex for the electric current one: reverberate, expand, outside, contract, and shrink. Our second improvement over PSO is to constitute a novel approach to diagnosing the convergence of the estimation. To exercise this, we propose three types of diagnostic measures: the local best-quartile method, the global best-variance method, and the local best-quartile-variance method.

The local best-quartile method uses the beginning and 3rd quartiles and the correlation structure of the population. Suppose θ^{one thousand} is the Due south × p matrix of the population (local best) of size S and the p parameters at kth iteration, i.e.,

$θ^{g} = [\begin{matrix} θ_{11}^{k} & \dots & θ_{1 p}^{k} \\ ⋮ & ⋱ & ⋮ \\ θ_{Due south i}^{k} & \dots & θ_{S p}^{k} \end{matrix}] due west i t h 10,$

where $θ_{i j}^{one thousand}$ is the local best of ith particle of jth parameter at kth iteration, 1 ≤ i ≤ S, 1 ≤ j ≤ p, 10 is the set of indices of each particle from ane to South and |X| = Due south.

We get-go calculate the commencement (lower) and 3rd (upper) quartiles, $Q_{1}^{one thousand j}$ and $Q_{3}^{k j}$ , for each parameter j at thouth iteration, where j = 1,two, ..., p, then obtain the reduced matrix ${\hat{θ}}^{k}$ using the offset and 3rd quartiles as follows:

${\hat{θ}}^{thousand} = [\begin{matrix} {\hat{θ}}_{11}^{thousand} & \dots & {\hat{θ}}_{1 p}^{thou} \\ ⋮ & ⋱ & ⋮ \\ {\hat{θ}}_{Thousand 1}^{thou} & \dots & {\hat{θ}}_{M p}^{m} \end{matrix}] w i t h 10^{k},$

where K = |X^yard | and Ten ^k ⊂X. In item, $X^{k} = \cap_{j = 1}^{p} {Ten}_{j}^{yard}$ and $X_{j}^{one thousand} = {i ∣ Q_{one}^{k j} \leq θ_{i j}^{grand} \leq Q_{three}^{k j}, i \in X}$ , j = 1,2, ..., p. And so the deviation between the showtime and third quartiles for each parameter is calculated based on ${\hat{θ}}^{k}$ , i.east.,

$d_{j}^{k} = ∣ Q_{1}^{k j} - Q_{three}^{k j} ∣, due west h due east r east j = 1, \dots, p,$

obtaining the maximum divergence of all the parameters as the following

In improver, the p × p correlation matrix of θ^k , i.eastward.,

is computed then its maximum and minimum eigenvalues, $λ_{yard a x}^{k}$ and $λ_{k i n}^{yard}$ , are estimated to calculate the ratio of ii eigenvalues, $ρ^{k} = ∣ \frac{λ_{m i due north}^{k}}{λ_{k a x}^{k}} ∣$ . If at to the lowest degree one parameter has $d_{j}^{m} = 0$ , then the eignvalues cannot exist obtained, and so we volition assign zero to ρ^thou in this case.

The global best-variance method considers the standard deviation of each judge of the parameters co-ordinate to the different window size. Suppose ψ^thousand is the g × p matrix consisting of the global all-time for each parameter up to kth iteration,

where $ψ_{j}^{i}$ is the global all-time of jth parameter at i the iteration and fifty^k is the vector of the loglikehood of each global best of size such equally 50^k = (l ¹, l ², ..., l^k ). Then the reduced matrix is obtained based on the user-divers window size, w. That is, the reduced matrix with west is

$ψ_{westward}^{k} = [\begin{matrix} ψ_{1}^{k - westward + ane} & \dots & ψ_{p}^{1000 - w + 1} \\ ⋮ & ⋱ & ⋮ \\ ψ_{one}^{k} & \dots & ψ_{p}^{one thousand} \end{matrix}],$

where k ≥ w < 0, and the reduced loglikehood vector is $l_{w}^{k} = (l^{m - due west + 1}, l^{k - west + 2}, \dots, l^{k})$ . Then nosotros compute the standard deviations for each parameter and loglikelihood:

$\begin{matrix} Southward D (ψ_{due west}^{k}) = \max_{j = ane, \dots, p} Due south D_{w} (ψ_{j}^{k}); \\ S D (l_{due west}^{thousand}) = \sqrt{V a r (l^{chiliad - w + i}, l^{grand - westward + two}, \dots, l^{chiliad}),} \end{matrix}$

where $Southward D_{westward} (ψ_{j}^{k}) = \sqrt{5 a r (ψ_{j}^{k - due west + 1}, ψ_{j}^{yard - westward + two}, \dots, ψ_{j}^{thou})}$ .

Besides, the local best-quartile-variance method is based on the standard deviation of the measures, d¹⁰⁰⁰ and ρ^{one thousand} , of the local all-time-quartile method. That is,

$\begin{matrix} South D (d_{w}^{k}) = \sqrt{5 a r (d^{k - west + 1}, d^{k - w + ii}, \dots, d^{k})}; \\ S D (ρ_{w}^{yard}) = \sqrt{V a r (ρ^{k - w + 1}, ρ^{k - w + 2}, \dots, ρ^{l})} . \end{matrix}$

Until now, nosotros innovate six unlike measures to diagnose the convergence of PSO at kthursday iteration: d^k, ρ^m , $S D (ψ_{w}^{k})$ , $Southward D (l_{w}^{chiliad})$ , $Southward D (d_{w}^{k})$ , and $S D (ρ_{w}^{thou})$ . Nosotros compare the performances of each arroyo in terms of the number of iterations to converge and the estimates at that iteration through the simulation in the next section.

The detailed procedure of the proposed LPSO is described as follows:

Stride 1. Initialization

For the parameter θ and the measurement mistake σ ^two, their populations are initialized randomly by their compatible distribution:

$x_{0} = (x_{0}^{1}, x_{0}^{2}, \dots, x_{0}^{s}, \dots, x_{0}^{S}) \sim U n i f (R (θ) \times R (σ^{2})),$

where southward = 1, ..., Due south; Due south is the size of population; $x_{0}^{s} = (θ_{0}^{south}, {(σ^{2})}_{0}^{s})$ is the sth particle for θ and σ ²; R(z) is the range of a random variable (vector) z; Unif stands for the multivariate uniform distribution.

Stride ii. A derivative-gratuitous local optimization-based estimation

For particle s at iteration k, the local optimization-based estimation, Nelder-Mead method, is rendered by an objective (loglikelihood) office. Let the current position of particle s at iteration k, $10_{1000}^{s} = (θ_{k}^{s}, {(σ^{2})}_{k}^{due south})$ exist an initial value for the estimation. Their estimates ${\hat{x}}_{grand}^{due south}$ are obtained and and then the current position $x_{k}^{s}$ is updated with the estimate ${\hat{x}}_{chiliad}^{s}$ which is converged and estimated by the local optimization algorithm, i.e.,

In example of failing to converge, there is no update on that particle.

Stride iii. Finding local and global all-time positions

The loglikelihood for particle s at iteration $k (l_{k}^{s})$ is calculated given the update current position $x_{chiliad}^{s}$ ,

$l_{yard}^{southward}$ is then compared to the best previous local and global best goodness of fits (i.e., $l_{l b e s t}^{s}$ and l_gbest ), and the current local and global bests are updated as follows:

Updating the local all-time position

${ten}_{l b e southward t}^{s} \leftarrow x_{k}^{due south} and l_{50 b due east s t}^{south} \leftarrow {fifty}_{thousand}^{s} i f l_{k}^{south} > l_{l b e s t}^{s},$

and there is no update otherwise;
Updating the global best position

$10_{one thousand b e due south t} \leftarrow x_{grand}^{s} and l_{g b e s t} \leftarrow l_{k}^{southward} i f {fifty}_{m}^{due south} > l_{thou b e south t},$

and there is no update otherwise.

Pace 4. Convergence

If iteration k reaches the pre-specified maximum Thou or all the particles are converged past the convergence diagnosis, the proposed algorithm LPSO will terminate. Otherwise it will go to the next step.

Footstep v. Updating the current positions

If the current iteration is not satisfied with one of the convergence criteria described, the current positions ${10_{k}^{s}}_{south = 1}^{S}$ will be updated to ${x_{k + 1}^{s}}_{s = ane}^{Southward}$ by (4) and (v). That is,

where $v_{chiliad + 1}^{due south}$ is divers in (4).

5. Simulation and Application

Simulation Experiment

The 2-compartment Four model described in Section 2 was employed for the simulation experiments. This model has five PK parameters θ = (logV1, logV2, logCL12, logVmax, logKm) and variance σ ^two for the measurement error. For the sake of simplicity, we fixed Dose and logV1 as one and zero, respectively, in society for the systemic concentrations to range naught to one. The identifiable and not-identifiable cases were used for comparing analysis. Every bit for the identifiable case, the true values for logVmax and logKm were zero and −ii.3, respectively, and, as for the non-identifiable instance, their values were nix and 15, respectively. We generated the simulated data without measurement error and its generated data are depicted in Figure 1. Equally a result, the parameters to judge are θ = (logVmax, logKm). In Figure i, we compare the nonlinear ODE based PK model with the linear ODE based PK model. We telephone call them linear and nonlinear in the sense that the intrinsic clearance is dependent to time. The linear ODE-based PK model uses the post-obit intrinsic clearance instead:

That is, the intrinsic clearance is the ratio of Vmax to Km and so that there is no influence from the concentration. In fact, the trace plots are much unlike from each other when the model is identifiable (blackness solid and dotted lines with circle), while both trace plots are very similar to each other in case of non-identifiable models (read solid and dotted lines with triangle) tin can be seen in Figure 1.

An external file that holds a picture, illustration, etc. Object name is nihms-533785-f0001.jpg

The trace plots between time and concentration for simulation experiments

The solid line is for the nonlinear ODEs and the dotted line is for linear ODEs. When log(Km) = –2.3, the open black circle is used and the red triangle is used when log(Km) = 15.

The constants of PSO were taken as (c1, cii, w_max, w_min, K) = (two,2,0.9,0.3,5000), and the number of particles of each parameter was 10 (for PSO) or 5 (for LPSO). The parameter boundaries are (–twenty, 20). The true values are θ^true = (0, –2.iii) for the identifiable case and (0,15) for non-identifiable instance. For both PSO and LPSO, the same seed number was used to generate the initial population. We besides utilise the five different cutoff values to consider PSO and LPSO as converging to a global optimum such as 0, 10^–v, x^–iii, 10^–ane, and 1. Afterward 5000 iterations, we can detect that both PSO and LPSO converge to the true value for the identifiable case, but the non-identifiable case doesn't converge to the truthful value for both PSO and LPSO equally shown in Table 1. In addition, we can run across that, in example of non-identifiable model, LPSO will not exist stopped if d^{one thousand} is used with a cutoff value less than or equal to one.

Table 1

The number of iterations to converge and its gauge for PSO and LPSO using the divergence of the first and third quartiles (d^m ) and the ratio of the eigenvalues (ρ^k) at kth iteration.

			PSO			LPSO			Truthful
		CV	Iter	${\hat{θ}}_{V m a x}$	${\hat{θ}}_{K grand}$	Iter	${\hat{θ}}_{V m a ten}$	${\hat{θ}}_{K m}$	θ _Vmax	θ _Km
Identifiable	*d^m*	0
		10^−five				57	0.0000	−2.3026
		10⁻ⁱⁱⁱ				9	0.0000	−2.3026
		10⁻¹				3	0.0000	−2.3026
		1	216	0.0000	−2.3026	3	0.0000	−2.3026

	*ρ ^k*	0							0	−two.3
		10^−five
		ten⁻ⁱⁱⁱ	293	0.0000	−two.3026	v	0.0000	−two.3026
		10⁻⁵	75	−0.0054	−2.3350	ii	0.0000	−2.3026
		ane	1	−0.6308	−16.9981	i	0.0010	−2.2970

		5000	0.0000	−ii.3026	5000	0.0000	−2.3026

Nonidentifiable	*d^m*	0
		x^−v
		ten⁻³
		x⁻¹
		1	59	1.8228	16.8226

	*ρ ^k*	0
		10⁻⁵				6	0.6985	15.6990	0	15
		10⁻³	98	ane.8005	sixteen.8005	5	−2.0000	13.0005
		10⁻¹	13	1.9787	xvi.9777	2	−2.6834	12.3145
		i	1	ane.1539	sixteen.0780	1	−nine.7187	v.2450

		5000	1.8948	sixteen.8948	5000	3.1095	eighteen.1095

Effigy 2 (a) and (c) display the trace plots for the global best of each parameter and their loglikelihood over iteration, while the trace plots for d^chiliad and ρ^k over iteration are depicted in Figure 2 (b) and (d), for identifiable and non-identifiable cases, respectively. We tin see that LPSO reaches steady land faster than PSO does in this figure. In Figure ane(c), the global best estimates of each parameter have bigger variation than its loglikelihood for both PSO and LPSO due to the non-identifiability. This is considering several estimates share the exactly aforementioned loglikelihood then that loglikelihood reaches steady state earlier for non-identifiable cases. In terms of the number of iterations to converge, ρ^k converges earlier than d^k , and LPSO converges faster than PSO, as shown in Table 1. As well, the estimates of each parameter are closer to the true values in example of LPSO. However, in instance of non-identifiable, although PSO and LPSO converge to certain estimates when ρ^thou or d^thou is less than a cutoff value, its estimate is far from the true value. Particularly, logVmax has much biased estimates for both PSO and LPSO in Table 1. Interestingly, LPSO converges to the estimates close to the true value already just after 1^st iteration in case of the identifiable model.

An external file that holds a picture, illustration, etc. Object name is nihms-533785-f0002.jpg

An external file that holds a picture, illustration, etc. Object name is nihms-533785-f0003.jpg

An external file that holds a picture, illustration, etc. Object name is nihms-533785-f0004.jpg

An external file that holds a picture, illustration, etc. Object name is nihms-533785-f0005.jpg

The trace plots between iteration and the global best and its log-likelihood

In each plot, the outset column is past PSO and the result of LPSO is the 2d column. The trace plots for the global best and its loglikelihood of logVmax and logKm at each iteration are depicted when (a) logKm = –2.3 and (c) logKm = fifteen where the dotted line indicates the true value. The trace plots for the difference of the start and third quartiles (|Q ₁ – Q ₃|) for logVmax and logKm and the ratio between the minimum and maximum eigenvalues ( $∣ \frac{λ_{m i n}^{grand}}{λ_{thou a ten}^{k}} ∣$ ) of correlation matrix at each iteration are depicted when (b) logKm = –2.3 and (d) logKm = 15.

Since PSO is a stochastic method and the direction of each particle is selected randomly, the loglikelihood is not increased monotonically unlike from a gradient based approach. In other words, although it is not converged yet, PSO volition stay in the same estimate or loglikelihood value, causing the users to consider PSO equally converging to a global optimum, peculiarly, if the showtime derivative (or the difference) of the previous and current estimates is employed. To overcome this upshot, we instead consider the standard deviation co-ordinate to the 5 different window sizes, w = x, xx, 30, xl, and fifty, as described in the previous section. Effigy 3 shows the trace plots for the standard deviation (SD) of the four diagnostic measures, $ψ_{westward}^{m}$ , $50_{westward}^{k}$ , $d_{w}^{m}$ , and $ρ_{w}^{k}$ . In the identifiable case, LPSO behaviors more stable than PSO in terms of the trace plots of each standard deviation, while PSO reaches steady country earlier than LPSO when the model is non-identifiable equally shown in Figure 3. This is because the management of each particle in LPSO is the same as that of a local optimum by its combined local optimization. If the model is non-identifiable, there are many estimates having the exactly same loglikelihood so that the local optimization volition finish to the dissimilar estimates according to the different initial values. For this reason, the trace plot of LPSO equally depicted in the right column of Figure iii(d) is non stable. However, PSO relies on the previous local and global bests so that it volition fluctuate less than LPSO. Tables 2 and three evidence the number of iterations to converge and the estimates according to the dissimilar window sizes (w) and cutoff values (CV). In case of the identifiable models, if the window size is large, it seems that all the methods tin can stop both PSO and LPSO when these algorithms converge to the true values for all the cutoff values. Still, if the window size becomes smaller, then PSO is frequently stopped before it converges to the truthful values, while LPSO e'er stops when information technology converge to the truthful values regardless of the diagnostic measures. This demonstrates that LPSO converges much faster than PSO. Overall, LPSO has a better performance than PSO in terms of the number of iterations to converge and the estimates, and $S D (ρ_{w}^{k})$ with the cutoff values less than x^–iii and the window size of ten seems to detect the correct iteration to stop both PSO and LPSO as can be seen in Tables ii and 3.

An external file that holds a picture, illustration, etc. Object name is nihms-533785-f0006.jpg

An external file that holds a picture, illustration, etc. Object name is nihms-533785-f0007.jpg

An external file that holds a picture, illustration, etc. Object name is nihms-533785-f0008.jpg

An external file that holds a picture, illustration, etc. Object name is nihms-533785-f0009.jpg

The standard departure plots for convergence

The standard deviation plots of the global best and its loglikelihood at each iteration are depicted when (a) logKm = –ii.3 and (c) logKm = 15 according to the five different window sizes (10, xx, 30, forty, l). The standard deviation plots for the departure of the first and 3rd quartiles (|Q ₁ – Q ₃|) and for the ratio between the minimum and maximum eigenvalues of correlation matrix at each iteration are depicted when (b) logKm = –2.3 and (d) logKm = 15. In each plot, the first cavalcade is by PSO and the consequence of LPSO is the 2d column.

Table two

The number of iterations to converge (Iter) and estimates for PSO and LPSO using the standard deviation of the global best ( $South D (ψ_{w}^{thousand})$ ) and of its loglikelihood ( $South D (l_{west}^{k})$ ) at kthursday iteration with different cutoff values (CV).

(a) PSO
	Window size (west)		ten			20			30			40			fifty
		CV	Iter	${\hat{θ}}_{V thou a x}$	${\hat{θ}}_{K m}$	Iter	${\hat{θ}}_{V m a x}$	${\hat{θ}}_{K k}$	Iter	${\hat{θ}}_{Five 1000 a x}$	${\hat{θ}}_{Thou thou}$	Iter	${\hat{θ}}_{V m a x}$	${\hat{θ}}_{K m}$	Iter	${\hat{θ}}_{V m a x}$	${\hat{θ}}_{Grand m}$
Identifiable	$Southward D (ψ_{w}^{k})$	0	29	−0.021	−2.434	78	−0.005	−2.335	981	0.000	−2.303	991	0.000	−2.303	1001	0.000	−2.303
		10⁻⁵	29	−0.021	−2.434	78	−0.005	−2.335	221	0.000	−2.303	230	0.000	−2.303	240	0.000	−2.303
		x⁻³	29	−0.021	−ii.434	78	−0.005	−2.335	144	0.000	−2.302	154	0.000	−2.303	163	0.000	−two.303
		10⁻¹	29	−0.021	−2.434	39	−0.017	−two.404	49	−0.007	−2.348	59	−0.005	−two.335	69	−0.005	−two.335
		1	17	−0.200	−xv.909	39	−0.017	−2.404	49	−0.007	−ii.348	59	−0.005	−2.335	69	−0.005	−2.335

	$Due south D (50_{w}^{m})$	0	29	−0.021	−2.434	78	−0.005	−2.335	580	0.000	−ii.303	665	0.000	−2.303	675	0.000	−ii.303
		10⁻⁵	eighteen	−0.200	−fifteen.909	39	−0.017	−two.404	49	−0.007	−ii.348	59	−0.005	−2.335	68	−0.005	−ii.335
		10⁻³	12	−0.200	−15.909	21	−0.021	−2.434	31	−0.021	−2.434	41	−0.007	−two.348	51	−0.007	−ii.348
		ten^−ane	10	−0.197	−sixteen.903	xx	−0.021	−ii.434	30	−0.021	−two.434	40	−0.017	−2.404	50	−0.007	−ii.348
		1	10	−0.197	−16.903	20	−0.021	−two.434	30	−0.021	−2.434	40	−0.017	−2.404	fifty	−0.007	−2.348

Nonidentifiable	$S D (ψ_{westward}^{yard})$	0	21	1.979	16.978	31	one.979	xvi.978	41	1.979	xvi.978	142	1.902	16.902	152	i.902	16.902
		10⁻⁵	21	i.979	16.978	31	1.979	16.978	41	1.979	16.978	142	ane.902	xvi.902	152	one.902	16.902
		10⁻³	21	one.979	sixteen.978	31	one.979	16.978	41	i.979	16.978	142	1.902	xvi.902	152	ane.902	16.902
		10^−one	20	1.979	16.978	29	1.979	16.978	38	1.979	16.978	47	one.931	16.930	55	i.931	16.930
		1	10	1.668	16.675	xx	1.979	16.978	30	ane.979	xvi.978	40	1.979	xvi.978	fifty	1.931	16.930

	$S D (l_{w}^{yard})$	0	21	1.979	xvi.978	31	1.979	xvi.978	41	i.979	16.978	142	1.902	sixteen.902	152	1.902	sixteen.902
		x^−v	ten	ane.668	16.675	20	1.979	xvi.978	30	one.979	16.978	xl	one.979	16.978	fifty	1.931	16.930
		10⁻³	10	i.668	16.675	20	one.979	16.978	xxx	one.979	sixteen.978	40	i.979	16.978	50	i.931	16.930
		10⁻¹	10	1.668	16.675	20	1.979	16.978	30	one.979	16.978	forty	1.979	16.978	l	1.931	16.930
		i	10	1.668	16.675	20	1.979	16.978	30	1.979	xvi.978	40	1.979	16.978	fifty	1.931	16.930

(b) LPSO
	Window size (westward)		10			20			30			twoscore			fifty
		CV	Iter	${\hat{θ}}_{5 m a x}$	${\hat{θ}}_{Km}$	Iter	${\hat{θ}}_{V m a x}$	${\hat{θ}}_{Km}$	Iter	${\hat{θ}}_{5 m a 10}$	${\hat{θ}}_{Km}$	Iter	${\hat{θ}}_{V yard a ten}$	${\hat{θ}}_{Km}$	Iter	${\hat{θ}}_{V m a x}$	${\hat{θ}}_{Km}$
Identifiable	$South D (ψ_{w}^{thou})$	0	31	0.000	−2.303	41	0.000	−2.303	51	0.000	−ii.303	106	0.000	−2.303	116	0.000	−2.303
		10⁻⁵	11	0.000	−2.303	21	0.000	−2.303	31	0.000	−two.303	41	0.000	−2.303	51	0.000	−2.303
		10⁻³	xi	0.000	−2.303	21	0.000	−2.303	31	0.000	−2.303	40	0.000	−2.303	50	0.000	−2.303
		x⁻¹	10	0.000	−ii.303	twenty	0.000	−ii.303	30	0.000	−ii.303	twoscore	0.000	−two.303	50	0.000	−ii.303
		1	10	0.000	−2.303	20	0.000	−ii.303	thirty	0.000	−2.303	40	0.000	−2.303	50	0.000	−2.303

	$Due south D ({fifty}_{w}^{k})$	0	31	0.000	−2.303	41	0.000	−ii.303	51	0.000	−2.303	106	0.000	−ii.303	116	0.000	−2.303
		x⁻⁵	10	0.000	−two.303	20	0.000	−2.303	30	0.000	−two.303	40	0.000	−2.303	l	0.000	−2.303
		10⁻³	10	0.000	−2.303	20	0.000	−2.303	30	0.000	−2.303	40	0.000	−2.303	fifty	0.000	−2.303
		10⁻¹	x	0.000	−two.303	20	0.000	−2.303	30	0.000	−2.303	40	0.000	−2.303	50	0.000	−2.303
		ane	x	0.000	−2.303		0.000	−ii.303	30	0.000	−2.303	40	0.000	−2.303	50	0.000	−2.303

Nonidentifiable	$S D (ψ_{due west}^{k})$	0	27	ii.778	17.778	37	ii.778	17.778	47	2.778	17.778	57	2.778	17.778	67	2.778	17.778
		10^−v	27	ii.778	17.778	37	2.778	17.778	47	2.778	17.778	57	two.778	17.778	67	2.778	17.778
		10⁻³	27	2.778	17.778	37	2.778	17.778	47	2.778	17.778	57	ii.778	17.778	67	2.778	17.778
		10⁻¹	27	2.778	17.778	37	ii.778	17.778	47	2.778	17.778	56	ii.778	17.778	66	2.778	17.778
		1	19	2.778	17.778	26	2.778	17.778	35	2.778	17.778	45	ii.778	17.778	54	ii.778	17.778

	$S D ({fifty}_{w}^{k})$	0	27	two.778	17.778	37	2.778	17.778	47	2.778	17.778	57	ii.778	17.778	67	2.778	17.778
		10^−v	10	4.271	19.271	20	2.778	17.778	30	2.778	17.778	twoscore	2.778	17.778	l	2.778	17.778
		ten⁻³	ten	4.271	19.271	twenty	two.778	17.778	thirty	2.778	17.778	40	two.778	17.778	50	2.778	17.778
		10⁻¹	ten	4.271	19.271	20	2.778	17.778	30	2.778	17.778	forty	2.778	17.778	50	ii.778	17.778
		1	10	4.271	19.271	twenty	2.778	17.778	30	2.778	17.778	40	2.778	17.778	50	2.778	17.778

Tabular array iii

The number of iterations to converge (Iter) and the estimates for PSO and LPSO using the standard divergence of the deviation of the outset and 3rd quartiles ( $Due south D (d_{w}^{k})$ ) and of the ratio of the eigenvalues ( $S D (ρ_{w}^{m})$ ) at kth iteration with dissimilar cutoff values (CV)

(a) PSO
	Window size (due west)		ten			twenty			thirty			40			50
		CV	Iterk	${\hat{θ}}_{V m a x}$	${\hat{θ}}_{Km}$	Iter	${\hat{θ}}_{V m a x}$	${\hat{θ}}_{Km}$	Iter	${\hat{θ}}_{Five m a 10}$	${\hat{θ}}_{Km}$	Iter	${\hat{θ}}_{V grand a x}$	${\hat{θ}}_{Km}$	Iter	${\hat{θ}}_{V chiliad a x}$	${\hat{θ}}_{Km}$
Identifiable	$South D (d_{westward}^{chiliad})$	0	106	−0.001	−2.307	207	0.000	−ii.303	245	0.000	−2.303	255	0.000	−2.303	265	0.000	−two.303
		x⁻⁵	106	−0.001	−ii.307	207	0.000	−2.303	245	0.000	−2.303	255	0.000	−2.303	265	0.000	−ii.303
		x⁻³	106	−0.001	−2.307	150	0.000	−2.302	214	0.000	−two.303	255	0.000	−2.303	265	0.000	−2.303
		10^−one	42	−0.007	−2.348	134	0.000	−2.303	143	0.000	−2.302	151	0.000	−two.302	160	0.000	−2.303
		i	19	−0.200	−15.909	xxx	−0.021	−2.434	39	−0.017	−2.404	114	0.001	−2.299	124	0.000	−2.303

	$Due south D (ρ_{west}^{m})$	0	1411	0.000	−2.303	1437	0.000	−2.303	1765	0.000	−2.303	1775	0.000	−2.303	1785	0.000	−2.303
		10^−v	318	0.000	−two.303	334	0.000	−2.303	347	0.000	−2.303	359	0.000	−two.303	369	0.000	−two.303
		x⁻ⁱⁱⁱ	106	−0.001	−two.307	146	0.000	−ii.302	156	0.000	−2.303	282	0.000	−two.303	292	0.000	−2.303
		10⁻¹	12	−0.200	−xv.909	41	−0.007	−2.348	51	−0.007	−2.348	64	−0.005	−2.335	75	−0.005	−two.335
		one	ten	−0.197	−sixteen.903	20	−0.021	−two.434	30	−0.021	−ii.434	xl	−0.017	−two.404	fifty	−0.007	−two.348

Nonidentifiable	$S D (d_{due west}^{k})$	0	92	1.842	xvi.842	151	ane.902	16.902	188	1.895	16.895	198	1.895	16.895	208	1.895	16.895
		ten⁻⁵	92	1.842	16.842	151	1.902	16.902	188	1.895	16.895	198	1.895	16.895	208	one.895	16.895
		10⁻³	92	1.842	16.842	151	1.902	16.902	188	one.895	xvi.895	198	1.895	16.895	208	1.895	16.895
		x⁻¹	31	ane.979	16.978	91	1.842	16.842	100	one.801	16.800	110	1.902	16.902	119	1.902	16.902
		1	21	1.979	16.978	31	1.979	16.978	42	1.931	16.930	61	1.823	16.823	75	1.872	16.872

	$South D (ρ_{due west}^{k})$	0	478	1.895	xvi.895	602	1.895	16.895	1154	1.895	sixteen.895	1164	one.895	xvi.895	1174	1.895	sixteen.895
		10⁻⁵	145	one.902	xvi.902	155	ane.902	16.902	196	1.895	xvi.895	206	i.895	xvi.895	216	1.895	sixteen.895
		10^−three	37	1.979	sixteen.978	47	1.931	16.930	62	1.823	16.823	72	1.872	16.872	82	1.872	sixteen.872
		10⁻ⁱ	21	i.979	sixteen.978	xxx	1.979	16.978	forty	1.979	16.978	50	1.931	sixteen.930	59	1.823	sixteen.823
		1	10	1.668	16.675	twenty	1.979	sixteen.978	30	1.979	16.978	twoscore	1.979	xvi.978	50	one.931	sixteen.930

(b) LPSO
	Window size (w)		10			20			thirty			40			fifty
		CV	Iter	${\hat{θ}}_{Five yard a x}$	${\hat{θ}}_{Km}$	Iter	${\hat{θ}}_{Five m a x}$	${\hat{θ}}_{Km}$	Iter	${\hat{θ}}_{V grand a x}$	${\hat{θ}}_{Km}$	Iter	${\hat{θ}}_{V m a ten}$	${\hat{θ}}_{Km}$	Iter	${\hat{θ}}_{5 m a x}$	${\hat{θ}}_{Km}$
Identifiable	$South D (d_{west}^{thousand})$	0	143	0.000	−ii.303	193	0.000	−2.303	203	0.000	−ii.303	213	0.000	−2.303	223	0.000	−2.303
		ten⁻⁵	62	0.000	−2.303	71	0.000	−2.303	81	0.000	−2.303	91	0.000	−2.303	101	0.000	−two.303
		x⁻³	12	0.000	−ii.303	22	0.000	−ii.303	32	0.000	−ii.303	42	0.000	−two.303	52	0.000	−2.303
		ten^−one	12	0.000	−2.303	22	0.000	−2.303	32	0.000	−2.303	42	0.000	−2.303	52	0.000	−2.303
		one	12	0.000	−2.303	22	0.000	−2.303	32	0.000	−2.303	42	0.000	−ii.303	52	0.000	−2.303

	$Southward D (ρ_{w}^{m})$	0	187	0.000	−ii.303	755	0.000	−2.303	812	0.000	−2.303	822	0.000	−2.303	832	0.000	−ii.303
		ten⁻⁵	187	0.000	−2.303	206	0.000	−2.303	209	0.000	−2.303	530	0.000	−2.303	540	0.000	−two.303
		x^−three	70	0.000	−2.303	94	0.000	−2.303	103	0.000	−2.303	112	0.000	−ii.303	116	0.000	−2.303
		10⁻¹	11	0.000	−2.303	21	0.000	−two.303	thirty	0.000	−ii.303	twoscore	0.000	−ii.303	50	0.000	−2.303
		i	10	0.000	−two.303	20	0.000	−2.303	30	0.000	−two.303	40	0.000	−2.303	50	0.000	−2.303

Nonidentifiable	$South D (d_{w}^{k})$	0	63	2.778	17.778	164	1.729	xvi.729	174	ane.729	sixteen.729	184	1.729	16.729	194	1.729	xvi.729
		10^−v	63	2.778	17.778	164	one.729	16.729	174	1.729	16.729	184	1.729	16.729	194	1.729	sixteen.729
		10⁻³	63	2.778	17.778	164	1.729	sixteen.729	174	one.729	16.729	184	ane.729	sixteen.729	194	1.729	16.729
		x^−ane	sixteen	2.202	17.202	164	1.729	xvi.729	174	i.729	sixteen.729	184	1.729	16.729	194	1.729	16.729
		1	fifteen	ii.202	17.202	49	2.778	17.778	60	2.778	17.778	108	one.729	xvi.729	149	ane.729	xvi.729

	$S D (ρ_{w}^{k})$	0	144	ane.729	16.729	173	1.729	16.729	183	1.729	xvi.729	326	2.211	17.211	336	two.211	17.211
		10⁻⁵	14	2.202	17.202	24	two.778	17.778	34	2.778	17.778	44	two.778	17.778	54	2.778	17.778
		10⁻³	14	ii.202	17.202	24	2.778	17.778	34	2.778	17.778	44	2.778	17.778	54	2.778	17.778
		10⁻¹	10	4.271	19.271	20	2.778	17.778	30	2.778	17.778	xl	2.778	17.778	50	2.778	17.778
		ane	ten	iv.271	xix.271	20	two.778	17.778	xxx	2.778	17.778	40	2.778	17.778	50	2.778	17.778

Effigy 4 shows the besprinkle-box plots for PSO and LPSO at 1^st, 500^th, and 5000^th iterations. As for the identifiable model, LPSO reaches closely the truthful value after 1^st iteration as depicted in the upper row of Figure 4 (a). Both PSO and LPSO display the bear witness of being non-identifiable after 500^th iteration in the non-identifiable instance as shown in the lesser row of Figure 4 (b). After 5000^thursday iteration, although PSO reaches the true value, it seems that logKm may still demand more iterations in the sense that there are a lot of outliers in its box plot in instance of the identifiable model as can be seen in the upper row of Effigy four(c). It is noteworthy that even though the pattern of the identifiable example of LPSO in Figure iv(c) is like to the non-identifiable, the ranges of x- and y-axis of each parameter are much narrower than those of the non-identifiable model. Upon investigating of these properties, the post-obit convergence diagnostics is proposed for LPSO.

An external file that holds a picture, illustration, etc. Object name is nihms-533785-f0010.jpg

An external file that holds a picture, illustration, etc. Object name is nihms-533785-f0011.jpg

An external file that holds a picture, illustration, etc. Object name is nihms-533785-f0012.jpg

The besprinkle-box plots betwixt logVamx and logKm for PSO and LPSO

The plots in the left and correct columns are for PSO and LPSO, respectively, and the first and second rows are for identifiable and non-identifiable cases. The solid lines in the plot bespeak the true values for each parameter.

Convergence diagnostics

If $S D (ρ_{w}^{g})$ is less than equal to the user-defined cutoff value (α) with the window size of w, LPSO will be considered as converged to a global optimum. Furthermore, if d^thousand is greater than the user-defined cutoff value (β), the model is considered as not-identifiable, where g is the number of iterations to converge which is identified past $S D (ρ_{due west}^{chiliad})$ . The general guideline for α and β is 0.001 and one, respectively.

Midazolam pharmacokinetic information analysis

We depict the analysis of a Midazolam (MDZ) pharmacokinetic data using the proposed approaches of convergence diagnostics with LPSO. MDZ is a benzodiazepine used to cause relaxation or sleep before surgery and to cake the memory of the procedure. It can be administrated in both oral and intravenous formulations. The MDZ PK study was conducted in the Full general Clinical Research Eye (GCRC) at Indiana Academy. Twenty-two subjects were recruited into this study. Blood samples for MDZ assays were nerveless in non-heparinized evacuated claret drove tubes at 0.five, 0.75, one, 1.v, 2, 4, 6, and 9 hours later on intravenously dosing MDZ (2.98 mg ~ 4.8 mg). We investigate iv of 24 subjects' clinical trial data for this study. Here the MDZ PK is causeless to follow a two-compartmental 4 model with the MM equation as described in Section ii. As well, this model has 6 log-transformed parameters equanimous of five PK parameters and variance, resulting in θ = (logV1, logVii, logCL12, logVmax, logKm, logσ ²). The constants of LPSO were set to (c ₁, c _ii, w_max, w_min, K) = (ii,ii,0.9,0.three,500) as the simulation studies did, and the number of particles of each parameter was fix to 3. The parameter purlieus for logσ ⁱⁱ is (–xx,0) and others are set to (–15, 15). The $South D (ρ_{w}^{m})$ is used every bit a stopping dominion with the cutoff value of α = 0.001 and the window size of w = 10, meaning that LPSO volition terminate if $Due south D (ρ_{10}^{one thousand}) \leq 0.001$ , where thousand is the current iteration.

The converged global optima for MDZ information are shown in Tabular array iv. The estimates of Km for Subject one and two are larger, while those for Subject three and four are close to zero, indicating that the parameter Km of these four MDZ information might exist non-identifiable. Namely, in instance of Bailiwick 1 and two, the approximate of Km is much larger than their concentrations, while the estimate of Km is much smaller than their concentrations. As a result, the estimates of Vmax for Subject 1 and ii become larger since only the ratio of Vmax to Km (Vmax/Km) is identifiable. The number of iterations to converge ranges from 67 to 388. Although the MDZ information set of Subject field three has the largest number of iterations to converge, it has the worst MSE (MSE = 3.8851), while Subject 1 has the smallest MSE (MSE = 0.0480) with the largest variance ( ${\hat{σ}}^{2} = 0.fifteen$ ). Note that MSE stands for the mean squared error between observed and predicted concentrations after log-transformed.

Table 4

The global optima for the parameters of MDZ information

The MSE stands for the mean squared error between observed and predicted concentrations after log-transformed.

Subject	V1	V2	CL12	Vmax	Km	σ ^two	MSE	The number of iterations
1	0.nineteen	0.28	0.xi	2593897	50278.35	0.fifteen	0.0480	133
2	5.01	19.55	eleven.45	406634.8	360.12	0.02	0.2455	191
3	4.77	51.07	27.87	36187.57	3.21E−7	0.02	3.8851	388
four	33.39	65.02	67.37	41628.6	iv.54E−vii	0.03	1.0600	67

Figure v displays several trace plots of the interpretation results of LPSO for four MDZ individual information from Field of study 1 to 4. For each subject field, the trace plots of global estimates of each parameter, the trace plots of loglikelihood and $South D (ρ_{10}^{k})$ , the trace plots of d^k , and the prediction plot are depicted in (a)-(d), respectively.

An external file that holds a picture, illustration, etc. Object name is nihms-533785-f0013.jpg

An external file that holds a picture, illustration, etc. Object name is nihms-533785-f0014.jpg

An external file that holds a picture, illustration, etc. Object name is nihms-533785-f0015.jpg

An external file that holds a picture, illustration, etc. Object name is nihms-533785-f0016.jpg

The results of estimation of LPSO for MDZ data

For each subject of four MDZ individual data, the trace plots of global estimates, the trace plots of the loglikelihood and the standard deviation of the ratio of eigenvalues ( $S D (ρ_{10}^{1000})$ ), the trace plots of the departure between Q1 and Q3 (d^k ), and the prediction plot are displayed in (a)-(d), respectively. The grey dotted horizontal line indicates the threshold of 1 to see whether a parameter is identifiable in (c).In (d), the MSE stands for the hateful squared fault between the log-transformed observation and the predicted concentration. The estimated Km value is indicated past the dotted grey horizontal line. The plots of Bailiwick i, 2, 3, and four are in the left-tiptop, the right-top, the left-bottom, and the right-bottom, respectively.

If we consider the trace plots of loglikelihood in Figure 5(a), all subjects reach the steady state afterwards ~10 iterations of LPSO, while $S D (ρ_{x}^{yard})$ notwithstanding fluctuates until information technology converges. On the other hand, compared to the trace plots of loglikelihood, the global optima for each private parameter still fluctuate as shown in Figure five(b). In Figure five(b), we can further observe several steady-land-like periods before reaching the convergence. For example, in case of Subject 1, all traces of the parameters take the steady-state period between 50 and 100 iterations before stopped. Therefore, if the stopping rule relies on the slope or deviation between 2 consecutive estimates, it is possible that the estimation volition stop between 50 and 100 iterations and and so give united states of america a local optimum instead of a global optimum. This demonstrates an advantage of the proposed convergence criteria over the gradient or difference based approaches.

The trace plots of d^k for logKm are larger than the cutoff of β = i across all the subjects, which indicate the non-identifiability every bit shown in Figure 5(c). In other words, according to the proposed Convergence Diagnostics, nosotros can make up one's mind that PK models of these four subjects are non-identifiable since there are i or more parameters having d^k greater than the cutoff of β = 1 until it converges. Nosotros can come across an interesting fact in Effigy 5(d). The MSEs of Subject 1 and 2 are smaller than those of Subject iii and four. As shown in Tabular array 4, the departure between Subject 1 and 2 and Bailiwick 3 and 4 is the size of estimates of Km. As for the showtime two subjects, their estimates of Km are large values, while the estimates of Km are near null for the last 2 subjects.

Overall, the clinical PK information analysis is consistent with the simulation studies and shows that the proposed LPSO and Convergence Diagnostics are able to not only diagnose the convergence of LPSO but also detect the identifiability.

6. Conclusion

A novel version of PSO is proposed with enhancing the convergence of the local best using a derivative-costless local optimization algorithm, which is called LPSO. In fact, the simulation studies and MDZ PK data analysis show that LPSO converges to a global optimum much faster than PSO does. Since PSO is a derivative-costless algorithm and a derivative-gratuitous local optimization is combined, the proposed LPSO becomes a derivative-free global optimization algorithm so that LPSO can be applied to the parameter estimation regardless of the identifiability. Furthermore, several convergence diagnostic measures are proposed and evaluated through both the simulation studies and clinical PK data analysis. Of these measures, using the maximum of the difference between the beginning and third quartiles and the standard difference of the ratio of the minimum and maximum of eigenvalues tin can detect when to stop LPSO besides as signal whether the model is identifiable or not.

Acknowledgements

Dr. Seongho Kim's research is partially supported by NSF Grant DMS-1312603 (Southward.Grand.). Dr. Lang Li's research is partially sponsored by NIH Grants R01 GMS 74217 and U54 CA113001 (Fifty.L.). The Biostatistics Core is supported, in part, by NIH Center Grant P30 CA022453 to the Karmanos Cancer Institute at Wayne State University.

Footnotes

Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early on version of the manuscript. The manuscript volition undergo copyediting, typesetting, and review of the resulting proof earlier it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the periodical pertain.

Conflict of Interest:

None declared

References

Beal SL, Sheiner LB. Estimating population kinetics. CRC Crit. Rev. Bioeng. 1982;viii:195–222. [PubMed] [Google Scholar]
Chang Grand. Monte Carlo simulation for the pharmaceutical manufacture: concepts, algorithms, and case studies. CRC Printing; Boca Raton, FL: 2010. [Google Scholar]
Eberhart R, Kennedy J. A New Optimizer Using Particle Swarm Theory.. Proceedings of the Sixth International Symposium on Micro Machine and Human Science; Nagoya, Nihon. 1995.pp. 39–43. [Google Scholar]
Engelbrecht AP. Computational Intelligence: An Introduction. John Wiley & Sons; 2007. [Google Scholar]
Godfrey KR, DiStefano JJ. Identifiability of model parameters in identifiability of parametric models. Pergamon Press; Oxford: 1987. [Google Scholar]
Godfrey KR, Fitch WR. The deterministic identifiablity of nonlinear pharmacokinetic models. J. Pharmacokin. Biopharm. 1984;12:177–191. [PubMed] [Google Scholar]
Hein GE, Niemann C. Steric class and specificity of a chymotrypsin catalyzed reactions, I. Journal of the American Chemical Society. 1962;84:4487–4494. [Google Scholar]
Hengl S, Kreutz C, Timmer J, Maiwald T. Data-based identifiability analysis of non-linear dynamical models. Bioinformatics. 2007;23:2612–2618. [PubMed] [Google Scholar]
Hidalgo ME, Ayesa Eastward. A numerical identifiability test for state-space models – application to optimal experimental design. Water Science & Technology. 2001;43:339–346. [PubMed] [Google Scholar]
Kennedy J, Eberhart R. Particle Swarm Optimization.. Proceedings of IEEE International Briefing on Neural Networks; Perth, Australia. Piscataway, NJ: IEEE Service Center; 1995. pp. 1942–1948. [Google Scholar]
Kim S, Li L. A novel global search algorithm on non-linear mixed-effect models using particle swarm optimization. Periodical of Pharmacokinetics and Pharmacodynamics. 2011;38:471–495. [PubMed] [Google Scholar]
Lecourtier Y, Lamnabhi-Lagarrigue F, Walter Eastward. In: Volterra and generating ability series approaches to identifiability testing, in identifiability of parametric models. Walter Due east, editor. Pergamon Printing; Oxford: 1987. [Google Scholar]
Lindstrom 1000, Bates D. Nonlinear mixed effects models for repeated measures data. Biometrics. 1990;46:673–687. [PubMed] [Google Scholar]
Metzler CM, Tong DDM. Computational problems of compartment models with Michaelis-Menten blazon elimination. J. Pharm. Sci. 1981;lxx:733–737. [PubMed] [Google Scholar]
Monod J. The growth of bacterial cultures. Ann Rev. Microb. 1949;3:371–394. [Google Scholar]
Nelder JA, Mead R. A simplex method for function minimization. Comput. J. 1965;7:308–313. [Google Scholar]
Morris MD. Factorial sampling plans for preliminary computational experiments. Technometrics. 1991;33:161–174. [Google Scholar]
Saccomani MP, Audoly S, D'Angio L. Parameter identifiability of nonlinear systems: the role of initial conditions. Automatica. 2003;39:619–632. [Google Scholar]
Shi YH, Eberhart RC. A Modified Particle Swarm Optimizer. IEEE International Conference on Evolutionary Ciphering; Anchorage, Alaska. May iv-nine, 1998. [Google Scholar]
Tong DDM, Metzler CM. Mathematical properties of compartment models with Michaelis-Menten blazon elimination. Mathematical Biosciences. 1980;48:293–306. [Google Scholar]
Vajda S, Godfrey KR, Rabitz H. Similarity transformation approach to identifiability analysis of nonlinear compartmental-models. Mathematical Biosciences. 1989;93:217–248. [PubMed] [Google Scholar]
Vonesh EF, Carter RL. Mixed-furnishings nonlinear regression for unbalanced repeated measures. Biometrics. 1992;44:673–683. [PubMed] [Google Scholar]
Wolfinger RD. Laplace's approximation for nonlinear mixed models. Biometrika. 1993;80:791–795. [Google Scholar]
Wynn HP, Parkin N. Sensitivity analysis and identifiability for differential equation models.. Proceedings of the 40th IEEE Conference on Decision and Control; Orlando, FL. 2001.pp. 3116–3121. [Google Scholar]
Yue H, Brownish M, He F, Jia J, Kell DB. Sensitivity analysis and robust experimental blueprint of a signal transduction pathway system. International Periodical of Chemical Kinetics. 2008;40:730–741. [Google Scholar]
Zhang T, Cai JD. A new chaotic PSO with dynamic inertia weight for economic acceleration problem. Sustainable Ability Generation and Supply, 2009. SUPERGEN '09. International Briefing on. :1–6. 6–7. [Google Scholar]

cookafteld.blogspot.com

Source: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3906649/