In drug discovery, efficacy and toxicity of the drug are the two of the most important endpoints in early phase trials. In phase I clinical trials, a main goal is to describe the dose-limiting toxicities ( D L T) and estimate the maximum tolerated dose ( M T D). Phase II clinical trials focus on ascertaining the efficacy of the drug and further evaluate the safety of the drug. For example, the most effective dose ( M E D ) that produces the maximum efficacy is a common quantity of interest. However, by itself it is not an adequate assessment of its usefulness without knowing its toxicity effects. For both economical and ethical and safety reasons, it is desirable to design a study that incorporates both efficacy and toxicity considerations at the onset.

To model trinomial responses above, [41] proposed the following constant slope proportional odds (PO) model for a patient’s response given dose x. log ( π 3 ( x ) / ( 1 − π 3 ( x ) ) ) = a 1 + b x log ( ( π 2 ( x ) + π 3 ( x ) ) / π 1 ( x ) ) = a 2 + b x , where π 1 ( x ), π 2 ( x ) and π 3 ( x ) correspond to, respectively, probabilities of observing “no reaction”, “efficacy without toxicity” and “toxicity” at dose x. For any dose x, ∑ i = 1 3 π i ( x ) = 1. The model assumes a constant effect of dose exists across the cumulative logits. Such assumption may be violated for a trinomial model, and is unlikely to be valid when there are more than three types of responses.

The continuation-ratio (CR) model is a more flexible alternative because it does not have the restrictive assumption as the PO model does as described in Chapter 9 in [2]. In the CR model, slopes of the two equations b 1 and b 2 can be constant or not. Assuming b 1 , b 2 > 0, the CR model is given by (1.1) log ( π 3 ( x ) / ( 1 − π 3 ( x ) ) ) = a 1 + b 1 x log ( π 2 ( x ) / π 1 ( x ) ) = a 2 + b 2 x . It is straightforward to show for this model, probability of “no reaction” is π 1 ( x ) = 1 ( 1 + e a 1 + b 1 x ) ( 1 + e a 2 + b 2 x ) ; probability of “efficacy without toxicity” is π 2 ( x ) = e a 2 + b 2 x ( 1 + e a 1 + b 1 x ) ( 1 + e a 2 + b 2 x ) ; and probability of “toxicity” is π 3 ( x ) = e a 1 + b 1 x 1 + e a 1 + b 1 x .

[16] found optimal approximate designs for estimating all parameters in the CR model and also selected function of the model parameters. Approximate designs are large sample designs or simply probability measures on the given design space. They are characterized by the number of points the design has, where the design points are and the proportion of observations to take at each of the point. The formal is called a locally D-optimal design and is found by maximizing the logarithm of the determinant of the information matrix. In contrast, a c-optimal design minimizes the asymptotic variance of the estimated function of interest obtained from the Delta’s method. In both cases, the optimization is over all designs on the given design space and the resulting design is termed locally optimal because the information matrix for a nonlinear model depends on the unknown model parameters. Hence to implement such designs, prior estimates or nominal values of the model parameters are required and they may come from a pilot study or from studies using a similar compound or drug. Clearly, poor estimates or incorrect nominal values for the model parameters can result in a different locally optimal design that is inefficient. To fix ideas, we focus on constructing locally optimal designs using metaheuristics but the technique can be easily amended to find other types of designs. The aim of this paper is to show that nature-inspired metaheuristic algorithms, which are already widely used in computer science and engineering, are both flexible and effective tools for finding multiple-objective optimal designs for clinical trials and tackling general medical problems. There are many appealing features of metaheuristics and the main ones are they virtually do not require technical assumptions for them to work well, they are generally fast, their codes are widely available and they are general-purpose optimization tools. Consequently, they are very flexible and can tackle very different types of optimization problems in increasingly many disciplines. Their meteoric rise in popularity in optimization, along with its reasons, are well documented; see for example, [49, 50].

Recently a class of algorithms called nature-inspired metaheuristic algorithms has proved very popular in the optimization literature. [49, 50] provided reasons for the rapid rise and interest in these algorithms. Early users of such algorithms to find exact optimal designs for linear models include [22, 4] who used an annealing algorithm to search for optimal designs for linear regression models, and [33] who used a genetic algorithm to construct exact D-optimal designs. Of particular note is the particle swarm optimization (PSO) introduced by [14] for tackling optimization problems. PSO is increasingly used across disciplines to solve hard optimization problems. PSO is essentially assumptions free and it searches in a simply and effective way. For example, unlike many algorithms, PSO does not require the objective function to be differentiable or convex and can solve non-convex high-dimensional optimization problems.

PSO is a metaheuristic optimization algorithm inspired from the way animals, such as birds and fishes, search for food. The birds fly continuously in the sky to look for food on the ground. Each has its own perception where the food is (local optimum) but it communicates with the rest and collectively decide as a flock where the food is (global optimum). Accordingly, each bird flies toward the global optimum in the direction of the local optimum (not giving up completely where it thinks the food is). Birds are referred as particles and each bird represents a candidate solution to the optimization problem. Velocities and locations of each bird are adjusted at each iteration and if and when the flock converges, the perceived global optimum is found. In order to efficiently identify the optimal points, we initiate a flock of birds in the pre-defined search space. Let X ( k ) be the locations of particles at the k-th iteration. Define t L ( k − 1 ) to be the locations with the best objective function values discovered by each particle before the k-th iteration, and t G ( k − 1 ) to be the locations of the best value found by the whole swarm before the k-th iteration. At the k-th iteration, the particles are updated by (2.1) X ( k ) = X ( k − 1 ) + V ( k ) , and (2.2) V ( k ) = w V ( k − 1 ) + c 1 R 1 ⊗ [ t L ( k − 1 ) − X ( k − 1 ) ] + c 2 R 2 ⊗ [ t G ( k − 1 ) − X ( k − 1 ) ] , where V ( k ) is the velocity of the particle. There are several parameters in (2.2). The inertia weight represents how active the birds are and is denoted by w. This parameter may be chosen to be a positive constant but more typically its value changes over iteration and eventually decrease to 0. The parameters c 1 and c 2 are two positive constants which are recommended to be 2, and R 1 and R 2 are two random vectors whose components are independently drawn from the uniform variate on [ 0 , 1 ]. In practice, the number of iterations and the swarm size are the most influential parameters in PSO. Large swarm size gives a better vision on the search area such that PSO could achieve the global optimum with a higher chance, and the more iterations allows the particles having search experience due to random perturbation. More details on PSO and the related metaheuristic optimization algorithms are available in [53]. Throughout, all computations were done on a PC with 3.50 GHz Intel(R) Core(TM) i7-4770K CPU.

The current literature on optimal designs is replete with single objective optimal designs. In practice, there are frequently two or more objectives of interest in a study and they may be of unequal interest to the researcher. There is some work on two and three-objective optimal design problems and we provide a brief review of them in this section.

When there are two objectives, the more important objective is the primary objective, and the other is the secondary objective. The sought dual-objective optimal design is one that optimizes the primary objective first, before optimizing the second objective and ensuring that the obtained design has a higher efficiency for the first objective. When there are multiple objectives, a similar idea applies; the objectives are first prioritized in terms of their importance and the sought design is the one that delivers higher efficiencies for the more important objectives. For example, in a model based dose response study, a researcher may wish to find a three-objective optimal design and in order of importance, the goals are to estimate the maximum tolerated dose, M T D, the most efficacious dose M E D, and the third goal is to maximize the precision for all model parameters’ estimates. The user-selected efficiencies for the sought design may be 90% and 80% for the first two objectives and subject to these constraints, the design should do as best possible under the tertiary objective. Does such a design exist and if it does, how do we find it?

For the simple case, when we have two-objective optimal designs and we have two linear and quadratic polynomials as competing regression models, [13] gave an analytical solution of such optimal designs. The dual-objective optimal design has to meet the user-specified efficiency requirement under the primary objective and is a constrained optimization problem. [13] showed that when the objectives are formulated in terms of concave functional of the information matrix, the constraint problem can be equivalently found under a simpler setup by optimizing a convex combination of the two (or more) concave function. This approach turns the dual-objective design problem into a single objective optimal design problem once the weights in the convex combination are fixed. Standard algorithms can then be applied to find the single-objective optimal design for each convex weight in the convex combination. The resulting designs are called compound optimal designs and Cook and Wong (1994) showed the class of compound optimal designs is the same as the class of constrained optimal designs for a dual-objective optimal design problem. To identify which compound optimal design coincides with the sought constrained optimal design with the user-specified efficiency under the primary criterion, they proposed plotting the efficiencies of each compound optimal design under the two criteria across all possible values of the weights from 0 to 1. The resulting plot is called an efficiency plot, which can then used to read off which convex combination weight will yield the sought constrained optimal design.

The simplicity of the above approach is valid for dual-objective optimal design problems. When there are multiple-objectives, the efficiency plot becomes high dimensional and it becomes difficult to visually identify the correspondence between which compound optimal design is the sought constrained optimal design. Unlike the two-objective optimal design problems, there is no explicit easy analytical description of the multi-objective designs, especially when the models are nonlinear and complicated. Our experience is that standard algorithms to find them also becomes problematic.

We found a variety of optimal designs for the flexible continuation-ratio (CR) model, which has great potential for dose finding studies because it simultaneously models both probabilities of observing efficacy and adverse effect without having to assume the correlation between them.

There is work in the literature that simultaneously study efficacy and toxicity by postulating bivariate parametric models. For example, to model dose-response curves, [24] proposed a bivariate logistic model which simultaneously measures toxicity (yes/no), and efficacy (yes/no), assuming these two outcomes are correlated. A similar idea was presented by [18] who proposed a dichotomized outcome pair ( y 1 , y 2 ) T from a bivariate normal model, where y 1 is the efficacy response and y 2 is the toxicity response. It follows that, there are four possible outcomes and the utility function of interest is π ( x ) = P r ( y 1 = 1 , y 2 = 0 | x ). This model is complicated and difficult to use in practice because it requires user to specify the correlation structure between toxicity and efficacy. Alternatively, we may combine the two outcomes corresponding to the presence of toxicity together and reduce the four types of responses into three. Doing so greatly simplifies the model by eliminating the correlation between toxicity and efficacy. In the resulting trinomial models, the responses of patients are now exclusively and exhaustively classified into three categories: “no effect” (neither toxicity or efficacy found); “efficacy” with no toxicity; and “adverse reaction” for toxicity regardless of whether efficacy is presented or not.

The optimal designs that we are interested to construct is a three-objective compound optimal design that provides efficient estimates for the most effective dose ( M E D), the maxim um tolerated dose ( M T D) for a user specified toxicity rate, and for all parameters in the CR model. We show that PSO can successfully find multi-objective compound optimal designs for the CR model. Further, we investigate the proper choice of weights for three optimal criteria in multi-objective designs under different parameters settings. By using efficiency plots, researchers and practitioners can construct the desired compound optimal design through appropriate weights combination of three optimal criteria in a more flexible and informative way.

[16] focused on finding the dose x that maximizes probability of “efficacy without toxicity” π 2 ( x ), which is selected as the most effective dose ( M E D). If b 1 = b 2 then M E D has a closed form given by M E D = − ( a 1 + a 2 ) / 2 b 1 , but if b 1 ≠ b 2 , there is no closed form solution of M E D. In their paper, we produce several locally c-optimal designs for estimating the M E D using different nominal parameter settings on a user-specified dose interval. The latter is likely more practically relevant than having an unrestricted dose interval. This feature is helpful because some algorithms work only on certain types of intervals only, especially when there are many interacting factors involved and the dose interval for each factor varies considerably citepxu.

In dose-response studies, side effects and drug toxicity can be very serious and need to be well controlled. One of the main targets in clinical trials is to find out the “dose that is closest to an acceptable level of toxicity”. The target dose is defined as Maximum Tolerated Dose ( M T D), which is the highest dose of a treatment agent that produces the dose limiting toxicity ( D L T) in a proportion ρ of patients. In the CR model, it follows that π 3 ( M T D , a 1 , b 1 ) = ρ . Here the probability ρ is user-specified and varies depending on what types of adverse effect is of interest. For life threatening side effects, ρ should be a small value. To find the M T D, practitioners usually start from a safe dose that produces a small ρ and gradually increase the dose to its highest acceptable level. The dose level at ρ is called lethal dose L D 100 ρ ([57]). For example, a common choice is the median lethal dose denoted by L D 50. For both economical and ethical reason, it is necessary to take both M E D and M T D (or L D 100 ρ ) into consideration at the design stage.

The worth of a design is measured of the quality of the Fisher’s information matrix. This matrix is proportional to the negative of the expectation of the second derivatives of the total log likelihood function with respect to the model parameters. This section derives the matrix for the CR model and makes clear the explicit nature of the objective functions.

The response from a subject assigned to dose x has one of three possible outcomes: y = ( y 1 , y 2 , y 3 ) T with ∑ i = 1 3 y i = 1. Its expected value is E ( y ) = π ( x ) T = ( π 1 ( x ) , π 2 ( x ) , π 3 ( x ) ) T as previously defined. Denote the parameters in the CR model by θ = ( a 1 , b 1 , a 2 , b 2 ) T and let δ x be the design that puts all its weight at dose x. To find the information matrix I ( δ x , θ ) of δ x for the CR model, we use the method developed by [58] for multi-nomial logistic models. Denote the left hand side of (1.1) by η ( x ) = log ( π 3 ( x ) / ( π 1 ( x ) + π 2 ( x ) ) ) log ( π 2 ( x ) / π 1 ( x ) ) log ( π 1 ( x ) + π 2 ( x ) + π 3 ( x ) ) and its right hand side by X θ = 1 x 0 0 0 0 1 x 0 0 0 0 a 1 b 1 a 2 b 2 . By the numerator-layout notation, we have ∂ η ( x ) ∂ π ( x ) = ∂ η ( x ) ∂ π 1 ( x ) ∂ η ( x ) ∂ π 2 ( x ) ∂ η ( x ) ∂ π 3 ( x ) = − ( π 1 ( x ) + π 2 ( x ) ) − 1 − ( π 1 ( x ) + π 2 ( x ) ) − 1 π 3 ( x ) − 1 − π 1 ( x ) − 1 π 2 ( x ) − 1 0 1 1 1 . Therefore the derivative of π ( x ) with respect to θ is derived by the chain rule G ( x ) = ∂ π ( x ) ∂ θ = ∂ π ( x ) ∂ η ( x ) ∂ η ( x ) ∂ θ = ( ∂ η ( x ) ∂ π ( x ) ) − 1 X = − π 1 ( x ) π 3 ( x ) − π 1 ( x ) π 2 ( x ) π 1 ( x ) + π 2 ( x ) π 1 ( x ) − π 2 ( x ) π 3 ( x ) π 1 ( x ) π 2 ( x ) π 1 ( x ) + π 2 ( x ) π 2 ( x ) ( π 1 ( x ) + π 2 ( x ) ) π 3 ( x ) 0 π 3 ( x ) 1 x 0 0 0 0 1 x 0 0 0 0

Since y T ∼ M u l t i n o m i a l ( 1 , π ( x ) T ), the likelihood function at dose x is L ( θ ; x ) = ∏ j = 1 3 π j y j ( x ) , and so the log-likelihood is l ( θ ; x ) = y T log π ( x ). By the chain rule, the score vector at dose x is ∂ l ( θ ; x ) ∂ θ = ∂ l ( θ ; x ) ∂ log π ( x ) ∂ log π ( x ) ∂ π ( x ) ∂ π ( x ) ∂ θ = G ( x ) T π 1 − 1 ( x ) 0 0 0 π 2 − 1 ( x ) 0 0 0 π 3 − 1 ( x ) y , and the unit information matrix I ( δ x , θ ) for the observation y at dose x is given by I ( δ x , θ ) = E ( ∂ l ( θ ; x ) ∂ θ ) ( ∂ l ( θ ; x ) ∂ θ T ) = G T ( x ) π 1 − 1 ( x ) 0 0 0 π 2 − 1 ( x ) 0 0 0 π 3 − 1 ( x ) (E y y T ) π 1 − 1 ( x ) 0 0 0 π 2 − 1 ( x ) 0 0 0 π 3 − 1 ( x ) G ( x ) . Since y T ∼ M u l t i n o m i a l ( 1 , π ( x ) T ), we have E y y T = E y E y T + V a r ( y ) = π 1 ( x ) 0 0 0 π 2 ( x ) 0 0 0 π 3 ( x ) , and so, I ( δ x , θ ) = G T ( x ) π 1 − 1 ( x ) 0 0 0 π 2 − 1 ( x ) 0 0 0 π 3 − 1 ( x ) G ( x ) = 1 x 0 0 0 0 1 x 0 0 0 0 T ( π 1 ( x ) + π 2 ( x ) ) π 3 ( x ) 0 0 0 π 1 ( x ) π 2 ( x ) π 1 ( x ) + π 2 ( x ) 0 0 0 1 1 x 0 0 0 0 1 x 0 0 0 0 = ( π 1 ( x ) + π 2 ( x ) ) π 3 ( x ) 1 x 0 0 x x 2 0 0 0 0 0 0 0 0 0 0 + π 1 ( x ) π 2 ( x ) π 1 ( x ) + π 2 ( x ) 0 0 0 0 0 0 0 0 0 0 1 x 0 0 x x 2 .

Locally c-optimal designs for estimating the M E D and M T D of the CR model with b 1 ≠ b 2 can be similarly derived. Recalling that M E D is the dose which maximizes π 2 the probability of having efficacy without toxicity. We derive M E D by solving d π 2 ( x ) d x = 0. After a few steps of simplification, we get M E D as the implicit solution of the following equation g ( x , θ ) = b 2 ( 1 + exp ( − a 1 − b 1 x ) ) − b 1 ( 1 + exp ( a 2 + b 2 x ) ) = 0 . To find the locally c-optimal design for estimating the M E D, we need to get its gradient ∇ M E D first. We follow [6]’s way to handle such equation g ( x , θ ) without an explicit solution for M E D as follows.

By the implicit function theorem, if the function g ( x , θ ) has continuous first derivatives, and M E D ( θ ) is continuous, then ∂ g ( x , θ ) ∂ θ | M E D = ∂ g ( x , θ ) ∂ x | M E D ∂ x ∂ θ | M E D . It follows that ∇ M E D T = ∂ x ∂ θ | M E D = ∂ g ( x , θ ) ∂ x | M E D − 1 ∂ g ( x , θ ) ∂ θ | M E D where ∂ g ( x , θ ) ∂ x | M E D = − b 1 b 2 ( exp ( − a 1 − b 1 M E D ) + exp ( a 2 + b 2 M E D ) ) , and ∂ g ( x , θ ) ∂ θ | M E D = − b 2 exp ( − a 1 − b 1 M E D ) − b 2 M E D exp ( − a 1 − b 1 M E D ) − ( 1 + exp ( a 2 + b 2 M E D ) ) − b 1 exp ( a 2 + b 2 M E D ) 1 + exp ( − a 1 − b 1 M E D ) − b 1 M E D exp ( a 2 + b 2 M E D ) T .

Thus, the locally c-optimal design for estimating the M E D maximizes − ∇ M E D T I − 1 ( ξ , θ ) ∇ M E D .

In the CR model, since log ( π 3 ( x ) / ( 1 − π 3 ( x ) ) ) = a 1 + b 1 x , a simple calculation shows M T D = 1 b 1 (logit ( ρ ) − a 1 ). Similarly, the locally c-optimal design for finding M T D is to maximize − ∇ M T D T I − 1 ( ξ , θ ) ∇ M T D , where ∇ M T D T = ∂ M T D ∂ θ = ( − 1 b 1 , a 1 − logit ( ρ ) b 1 2 , 0 , 0 ). In the traditional dose finding designs, M T D is usually estimated using a 3+3 design or set the MTD level as the dose level at which two of six patients experienced toxicity, see [20] for example. Therefore in what is to follow, we set ρ = 0.3, but other values can be directly used as well.

In what is to follow, we treat a design ξ as an approximate design, i.e. a probability measure on a given dose space. If N is the total number of subjects in the study and ξ has support at x 1 , … , x k with corresponding weights w 1 , … , w k and w 1 + ⋯ + w k = 1, then the implemented design is to take N w i at x i , i = 1 , … , k after rounding each N w i to a positive integer [ N w i ] and subject to [ N w 1 ] + ⋯ + [ N w k ] = N. In particular, this means that for a K-point approximate design with weights w i at dose x i , i = 1 , 2 , … , K, the normalized information matrix for the design ξ is I ( ξ , θ ) = ∑ i = 1 K w i I ( δ x i , θ )

[26] proposed working with approximate designs after noting that if the criteria are all formulated as convex or concave functional of the design information matrix, they are easier to construct, study and check for their optimality. The latter is accomplished via an equivalence theorem for each convex or concave functional. They all have the form but are different for different criteria. Recalling the directional derivative of Ψ ( I ( ξ , θ ) ) = log | I ( ξ , θ ) | at ξ in the direction of δ x is ψ ( x , ξ , θ ) = trace [ I ( δ x , θ ) I − 1 ( ξ , θ ) ] − P and P is the number of parameters in the mean function. For example, if b 1 ≠ b 2 in the CR model, P = 4. Applying standard convex analysis results for the D-optimality criterion, we conclude the following statements are equivalent ([26]):

1. The design ξ ∗ maximizes Ψ ( I ( ξ , θ ) );

2. The design ξ ∗ minimizes max x ∈ χ ψ ( x , ξ , θ );

3. max x ∈ χ ψ ( x , ξ ∗ , θ ) = 0, and it achieves its maximum at the support points of the design.

This is the equivalence theorem for D-optimality that allows us to use (3) to check optimality of any design ξ in practice. Each convex or concave functional has a unique equivalence theorem but they all have the same form.

In the same paper, [26] also showed the equivalence theorem to verify the optimality of a locally c-optimal designs for estimating M E D or M T D. Both can be subsumed under the criterion Ψ ( I ( ξ , θ ) ) = trace [ A I − 1 ( ξ , θ ) ] when A is chosen appropriately. For estimating the former, the directional derivative of Ψ ( I ( ξ , θ ) ) = − ∇ M E D T I − 1 ( ξ , θ ) ∇ M E D = -trace [ ∇ M E D ∇ M E D T I − 1 ( ξ , θ ) ] at ξ in the direction of δ x is ψ ( x , ξ , θ ) = trace [ I ( δ x , θ ) I − 1 ( ξ , θ ) ( ∇ M E D ∇ M E D T ) I − 1 ( ξ , θ ) ] − trace [ ( ∇ M E D ∇ M E D T ) I − 1 ( ξ , θ ) ] . When ξ is optimal, ψ ( x , ξ , θ ) is bounded above by 0 with equality at the support points. Similarly we can derive the equivalence theorem for locally c-optimal design to estimate other quantities such as M T D. Because equivalence theorems confirm the optimality or non-optimality of the design, these theorems are generally referred to also as checking conditions and the function ψ ( x , ξ , θ ) is referred as the sensitivity function of the design ξ.

In [16], only M E D was of interest and the authors focused on finding locally c-optimal design for estimating the M E D. This study extends their work to construct a design that efficiently estimates M E D and M T D simultaneously. Other parameters such as b 1 and b 2 determine the slopes of three possible curves and so it is important to make sure the generated optimal design provides good estimates for all parameters in the CR model. Accordingly, in addition to estimating M E D and M T D, we incorporate D-optimality as well. For a given λ, we construct a three-objective compound optimal design by maximizing the weighted sum of three concave functionals Ψ ( ξ | λ ) = λ 1 Ψ 1 ( I ( ξ , θ ) ) + λ 2 Ψ 2 ( I ( ξ , θ ) ) + λ 2 Ψ 3 ( I ( ξ , θ ) ) where λ = ( λ 1 , λ 2 , λ 3 ) T , λ i ’s ∈ [ 0 , 1 ] and ∑ i = 1 3 λ i = 1. Here Ψ 1 ( I ( ξ , θ ) ) is the log c-efficiency of the design ξ for estimating the M T D Ψ 1 ( I ( ξ , θ ) ) = log ( ∇ M T D T I − 1 ( ξ M T D , θ ) ∇ M T D ∇ M T D T I − 1 ( ξ , θ ) ∇ M T D ) = log ( e 1 ) , Ψ 2 ( I ( ξ , θ ) ) is log c-efficiency of the design ξ for estimating the M E D Ψ 2 ( I ( ξ , θ ) ) = log ( ∇ M E D T I − 1 ( ξ M E D , θ ) ∇ M E D ∇ M E D T I − 1 ( ξ , θ ) ∇ M E D ) = log ( e 2 ) , and Ψ 3 ( I ( ξ , θ ) ) is log D-efficiency of the design ξ for estimating all parameters in the model Ψ 3 ( I ( ξ , θ ) ) = log ( | I ( ξ , θ ) | | I ( ξ D , θ ) | ) 1 P = log ( e 3 ) .

The checking condition for a three-objective compound optimal design is just a weighted sum of the checking conditions for the optimal design for each objective, i.e. we have ξ is optimal for the three-objective problem if and only if for all x ∈ χ λ 1 trace [ I ( δ x , θ ) I − 1 ( ξ , θ ) ( ∇ M T D ∇ M T D T ) I − 1 ( ξ , θ ) ] trace [ ( ∇ M T D ∇ M T D T ) I − 1 ( ξ , θ ) ] + λ 2 trace [ I ( δ x , θ ) I − 1 ( ξ , θ ) ( ∇ M E D ∇ M E D T ) I − 1 ( ξ , θ ) ] trace [ ( ∇ M E D ∇ M E D T ) I − 1 ( ξ , θ ) ] + λ 3 trace [ I ( δ x , θ ) I − 1 ( ξ , θ ) ] P ≤ 1 , with equality at the support points.

We apply standard PSO to find the locally c-optimal design for estimating the L D 30 in the CR model with nominal values a 1 = − 3.3 and b 1 = 0.5. Here and elsewhere, we use 100 particles and 1000 iterations to find the locally optimal design. The PSO-generated design is supported at a single dose at 4.905, which can be directly verified by checking L D 30 = 1 0.5 (logit ( 0.3 ) + 3.3 ) = 4.9054. Similarly, we use PSO to find different types of locally compound optimal designs using various nominal values for the CR model and note that all generated designs have four or fewer points. The three sets of nominal values we have are selected to generate mean response curves that we expect, i.e. the maximum efficacy varies from 0.3 to 0.9, see Figures 1(a), (b) and (c). In each figure, probability curves are shown in the left plot, and the middle one is the equivalence plot for compound optimal designs in the unrestricted dose space. Table 1 shows selected compound optimal designs with equal weights λ i ’s=1/3 for all 3 criteria on the unrestricted designs space. We observe that generated optimal designs are supported at three or four points, and the weights at these support points are unequal.

Invariance of the optimal designs to linear transformation of the design is a desirable property for an optimal design. D-optimal designs for linear models are examples of optimal designs with such a property, but A-optimal designs, that minimize the average variance of all the estimated parameters, do not have such a property, i.e. an A-optimal design on a new interval cannot be deduced from the A-optimal design on the origin interval. Interestingly, some algorithm may work on a design space but not on another design space, especially when the latter has are many interacting factors and each factor is defined on intervals with very different widths; see [52] for details and examples when they found locally D-optimal designs for a logistic model with many factors and all pairwise interactions are included.

To further investigate the influence of dose space on the compound optimal design, we verified the algorithm was still able to find compound optimal designs on an unrestricted space that may not be symmetric about 0 for the 3 sets of nominal values for the model parameters. For example, Table 2 displays the PSO-generated compound optimal designs on the asymmetric dose interval [−2,7] for all the three cases. Interestingly, the previous four point optimal designs become three points, with one or two points located on the boundaries of the restricted design space. We also observe that the middle support point is not the same as the ones in optimal designs on unrestricted dose space. The plots on the right hand side of Figures 1(a), (b) and (c) display the corresponding sensitivity functions and visually confirm that all the PSO-generated three-objective designs are optimal among all designs on [−2,7].

In our PSO implementation, for each fixed λ in [0,1], we use the default parameter settings in the standard PSO to find the three-objective compound optimal design ξ λ for the CR model. Particles that wander outside of the dose space are pulled back to its nearest boundary [ l b , u b ]: z i ( t ) = u b if z i , ( t ) > u b l b if z i ( t ) < l b z i ( t ) otherwise where z i ( t ) is the particle i’s position at time t. The rationale for our strategy is that optimal designs frequently have support points at the boundary of the dose space, see [29] for example. Results from a small numerical experiment not included here suggest that generally, PSO with boundary repair outperforms PSO with random repair in searching for locally D- and c-optimal designs. Using equivalence theorems, we were also able to verify that PSO-generated locally designs under the D-optimal and c-optimality criteria for the CR model were optimal for the different sets of nominal values in [16].

Can we find a three-objective optimal design ξ for the CR model such that its D-efficiency e 3 ( ξ ) is maximized subject to both its c-efficiencies e 1 ( ξ ), e 2 ( ξ ) for M T D and M E D are equal to or greater than 0.9? Such a design may or may not exist depending on how competitive the criteria are and also the efficiency requirements in the constraints. If such a design exists, what values for the λ i ’s should be chosen in the compound optimal design?

To answer the question we need to investigate the impact of different combinations of λ 1 and λ 2 on the efficiencies of the compound optimal design ξ λ relative to the three criteria. We first discretize λ 1 and λ 2 using a grid size of 0.05, and find the corresponding compound optimal design ξ λ for each pair of λ 1 and λ 2 combination subject to λ 1 + λ 2 ≤ 1. The efficiencies e i ( ξ λ ) under the ith criterion are displayed in Figure 2, 3 and 4 for different nominal values sets. In these figures, the upper left plot shows the c-efficiencies of e 1 ( ξ λ ) as measured by the locally c-optimal design ξ M T D for M T D; the upper right plot shows c-efficiencies of e 2 ( ξ λ ) relative to the locally c-optimal design ξ M E D for M E D; the bottom left plot shows D-efficiencies of e 3 ( ξ λ ); bottom right shows minimal efficiencies of all three efficiencies.

These efficiency plots provide an answer to the question posed at the beginning of the section. The answer is no, at least for the three sets of nominal values we show here; there is no compound optimal design ξ λ that maximizes D-optimality with subject to e 1 ( ξ λ ) ≥ 0.9 and e 2 ( ξ λ ) ≥ 0.9 no matter what combination of λ i ’s we choose. By observing the efficiency plots of the three criteria under different nominal values sets, we find e 1 ( ξ λ ) and e 2 ( ξ λ ) seem to be very competitive with each other; see Figures 2, 3 and 4). As expected, when one criterion gains in value, the other loses, and vice versa. Competitive criteria will result in a substantial loss for a small gain in the other, and vice versa.

For each set of the nominal values, Table 3 displays the combination of λ i ’s that produces the maximum of the minimum efficiencies. For example, λ 1 = 0.55, λ 2 = 0.35, λ 3 = 0.1 is the weights combination that produces the maximum of the minimum efficiency as 0.63 for all three criteria. That means with this choice for λ, we can find a compound optimal design ξ λ that maximizes the precision for estimating all parameters and at the same time has at least 63% efficiency for estimating the M E D and M T D.

In addition, we observe that for all the three nominal values sets, the compound optimal designs maintain high D-efficiency for most of the combinations of λ i ’s. This indicates for those nominal values that consistently produce high D-efficiency for different criteria weight λ i ′ s, we may reduce the three-objective compound optimal design to two-objective for estimating the M T D and M E D only so that the optimal design problem is simplified.

We close this section by noting that the above approach of finding locally optimal designs can be directly extended to finding Bayesian optimal designs with one or multiple objectives. Other methods for finding Bayesian optimal designs for Phase I trials and do not require that each objective be a convex function of the information matrix are available. Some seminal work and notable work in this direction are [42, 40, 56], among many others.

As an application, we revisit the clinical trial in [56] to identify the biologically optimal dose for a biologic agent to boost haemoglobin (HbL) levels for patients whose HbL were below the normal range. After the drug administration, the outcomes of patients could be categorized as ‘no response’, implying the drug had no effect, ‘success’, implying that the HbL levels were back to the normal range or ‘toxicity’, meaning that patients were over-stimulated and their HbL were raised beyond the normal range. The probabilities of these event sum to unity and the authors also chose the continuation ration model because of its flexibility and particular usefulness for modelling the three probabilities realistically. The probabilities have unknown parameters and so the information matrix of a design also has unknown parameters in them. Consequently, the design criterion which is formulated as a convex function of the information matrix also contains unknown parameters. Like, [42, 40], [56] also used a Bayesian approach and so they specified prior densities for the unknown parameters and averaged them out in the design criterion before finding a design to optimize the criterion. Locally optimal designs used in this paper can be viewed the same except that the prior densities are all one point priors. Thus depending on the situation, the practitioner may prefer a Bayesian optimal design or a locally optimal design. In the latter case, a single best guess for each parameter is needed and the prior elicitation process may be easier for the practitioner to do so. Frequently best guesses come from prior experience with a similar drug or from a pilot study to estimate the model parameters; see [10] for implementation details of a locally optimal design.

When tackling design or any optimization problems, especially high dimensional ones, where the solution is unknown, it is always not clear which metaheuristic algorithm to use. Solely relying on results from a single metaheuristic algorithms is risky since they do not guarantee convergence to the global optimum. It is thus important to run several metaheuristic algorithms and observe whether they all produce similar solution. If it is the case, the confidence in our solution is increased.

We now briefly compare performance of PSO with two other algorithms. We choose to compare CPU time required by PSO to find optimal designs with two of its competitors: Differential evolution (DE) and Cocktail algorithm (CA). Some of the models we used in the comparisons are listed in Tables 4 and 5. The compartment model is taken from [5] and it has 3 parameters and some of the interesting features to estimate were area under the curve (AUC), the time to maximum ( t m a x ) and the model parameters. The Hill model is the 4-parameter logistic model and the bivariable linear model has two additive factors. All have normally and independently distributed errors with constant variance. In the following tables, the numbers in parentheses are flocks size and the maximum number of iterations allowed for PSO and the single number in parentheses represents the number of points in the discretized dose interval

DE is a metaheuristic algorithm like PSO, which is a general purpose optimization algorithm. We use it to find locally D-optimal design for models ( a ) to ( e ) in Table 4, and locally c-optimal designs for estimating the three quantities ( f ) to ( h ) in Table 5. In this section, we worked with known optimal designs and the generated designs by PSO, DE and CA are considered to be optimal when all their design points and their weights are the same as the known optimal designs rounded to the 4 decimal places.

CA is a deterministic algorithm, which means that the same input results in the same answer, which in our case, results in the same optimal design. CA requires the design space to be discretized into a finite number of dose levels to approximate the continuous space. This approach can limit the ability of CA to find the optimal design because the discretized points may not correspond to a design point of the optimal design. However, our experience is that the CA-generated designs can usually still be highly D-efficient (with D-efficiency greater than 90 %). In addition, we observe from Table 6 that the number of support points in the PSO-generated designs tend to be smaller that those generated by CA. Thus the simpler PSO-generated designs are attractive and also likely less costly to implement because of the fewer design points.

Cocktail algorithm (CA) is a deterministic algorithm which always produces the same output given one input. For PSO, we still use the same tuning parameters in section 5.2.1 such that the algorithm produces the locally D-optimal at 90% chance. Our empirical experience with CA is that although its generated designs can be very efficient even when the grid size n is as small as 10, they may have many more support points than the true optimal design does. For example, we compare the locally D-optimal designs for bivariable linear model d ) generated by PSO and CA in Table 6. The numbers in the parenthesis after PSO are flock size and maximum iteration number; the numbers in the parenthesis after CA are the grid size used in different optimization problems. PSO is able to find the true optimal design supported at 6 points with D-efficiency 1. Both of the CA generated designs have D-efficiency greater than 0.99, but they have two more support points near ( 0 , 0 ) T and ( 0 , 1 ) T due to the discretized dose space. In biomedical researches such as a dose response study, practitioners usually prefer a simple design with fewer support points to a complicated design with a lot of points on the condition that both designs provide almost the same efficiency. Therefore in addition to the average CPU time,we use the number of support points in the generated design as another comparison criterion to measure the accuracy of CA generated design.