In this section, the general case of an optimal dosing task regarding efficacy and safety targets is presented and formulated as a state-constrained OCP. Methods for solving such are discussed.

Typically, a PMX model is a system of parameter- and dose-dependent ordinary differential equations,

$$\begin \begin \frac \,y(t)&=f(t, y(t), \theta , ), \qquad y(t_0) = ( \theta ), \end \end$$

(5)

defined in the time interval \( t \in [t_0,t_f]\) where \(f :[t_0,t_f] \times \mathbb ^n \times \Theta \times \mathbb ^m \rightarrow \mathbb ^n\) describes the mechanism and y the model state. Further, \(\theta \) denotes the model parameter in the set \(\Theta \) of admissible model parameters. The doses \(=( _1, \ldots , _m ) \in \mathbb ^m\) are administered at fixed dosing time points according to a fixed dosing scenario. For given \(\), state equation Eq. (5) admits a unique solution \(y() :[t_0,t_f] \rightarrow \mathbb ^n\) (under smoothness assumptions on the right-hand side, see [1]).

Please note that Eq. (5) is a slightly simplified notation compared to [2], as we do not distinguish between PMX mechanism and inflow resulting from drug administration anymore, both are now contained in f in Eq. (5). Further, model parameter \(\theta \) needs to be reasonably estimated, characterizing an individual or typical representative of a population, and is fixed during the optimization.

The aim is to compute optimal doses regarding efficacy and safety targets. Mathematically, this means minimizing a so-called cost function \(J :\mathbb ^m \rightarrow \mathbb \) subject to a state constraint,

$$\begin g (y(t,)) \le 0, \qquad \hbox t \in [t_1, t_2], \end$$

(6)

with a function \(g :\mathbb ^n \rightarrow \mathbb \), where \(t_1\) indicates the start and \(t_2\) the end of the time interval where the state constraint holds. Note that \(g(y(t,)) = 1-N(t,)\) recovers Ineq. (1) with \(t_1 = 12, t_2 = 32\). Doses \(\) are called feasible, if associated model state \(y()\) is feasible, i.e., if it satisfies the state constraint.

The state-constrained OCP reads

$$\begin \min J() \hbox \, y() \hbox (5), \\ \, y() \hbox (6), \\ \, \hbox _ \le \le _, \end\right. } \end$$

(7)

where \(_, _ \in \mathbb ^m\) are a lower and upper bound on the doses, respectively.

Typically, both cost function and state constraint reflect the therapeutic goal given by efficacy and safety targets. In the motivational example, the cost function describes an efficacy and the state constraint a safety target. In other examples, as we will see in the “Results” section, this can be vice versa. The OCP structure in Eq. (7) allows for a top priority target (i.e., safety or efficacy) characterized by a state constraint, whereas the secondary target, described by an appropriate cost function, is reached as good as possible without violating the state constraint.

An approach to solve problem Eq. (7) are penalty methods [6, 7]. A state-constrained OCP is transformed into an unconstrained OCP depending on a penalty parameter. The objective function of this unconstrained OCP is given as sum of the cost function J of the state-constrained OCP and an additional penalty function P. A penalty function measures the violation of the state constraint Ineq. (6) multiplied by a penalty parameter \(\rho > 0\), e.g.,

$$\begin P(,\rho ) = \frac \int _^ \max \))\}^2 \; dt. \end$$

(8)

The resulting unconstrained OCP then reads

$$\begin \begin \qquad&\min J() + P(,\rho ) \\&\hbox \, y() \hbox \\ \, \hbox _ \le \le _. \end\right. } \end \end$$

(9)

The idea of penalty methods is to solve a series of unconstrained OCPs Eq. (9) with increasing penalty parameter \(\rho \), whose solutions converge towards the solution of the original state-constrained OCP Eq. (7) as \(\rho \rightarrow \infty \) [7]. In the optimization in NONMEM, we fix \(\rho \) to a large value.

To calculate cost function value (CFV) and penalty function value (PFV), we introduce additional state variables \(y_\) and \(y_\), in the same way as when calculating the AUC of a drug from its concentration over time. In the motivational example, we have

$$\begin \frac \,y_(t)&= W(t,) \; \hbox t \in [t_0,t_f], \qquad y_(t_0) = 0, \end$$

which yields the CFV

$$\begin J() = y_(t_f) = \int _^ W(t,)\; dt, \end$$

and analogously

$$\begin&\frac \,y_(t) = \frac \max \)\}^2, & t \in [t_1,t_2] \\ 0, & \hbox \end\right. } \\ &\, \quad y_(t_0) = 0 \end$$

providing the PFV

$$\begin P(,\rho ) = y_(t_f) = \frac \int _^ \max \)\}^2 \; dt. \end$$

Several state constraints can be included in the penalty method utilizing several penalty functions. Moreover, state constraints of another type may appear, e.g., some may act only at a specific time point instead of an interval, i.e., \(t_1 = t_2\) in Ineq. (6). Then, the integral in Eq. (8) reduces to evaluation of the integrand at that specific time point. In the motivational example, an additional state constraint could be to demand the neutrophil level at final time \(t_f\) above the threshold of 3, expressed by the penalty function

$$\begin P_2(D,\rho _2) = \frac \max \)\}^2 \end$$

(10)

and the unconstrained OCP

$$\begin \begin \qquad&\min J() + P(,\rho ) + P_2(,\rho _2) \\&\hbox \, y() \hbox \\ \, \hbox _ \le \le _. \end\right. } \end \end$$

Another example for a different disease could be, e.g., to require a trough concentration below a certain threshold as a safety target, or a maximal efficacious concentration above a certain threshold as an efficacy target.

While the concept of penalty methods allows for several state constraints, each included with their own penalty function and penalty parameter, this increases the complexity of the problem. It might occur that the state-constrained OCP has no feasible solution, i.e., no model state satisfying all state constraints exists. In this case, the underlying efficacy and safety targets cannot be realized simultaneously.

In all PMX optimal dosing examples studied so far, it seemed sufficient to apply the proposed penalty method. However, a more advanced approach for solving state-constrained OCPs Eq. (7) are safeguarded augmented Lagrangian methods [7,8,9], which will now be briefly discussed. An additional variable \(\lambda :[t_1,t_2] \rightarrow \mathbb \) is introduced which approximates the Lagrange multiplier associated with the state constraint.

A popular penalty function in augmented Lagrangian methods is the Powell–Hestenes–Rockefellar [8, 10] penalty function

$$\begin \begin P_(,\lambda ,\rho )= \frac \int _^&\max \))\}^2 \\ &- \lambda (t)^2 \; dt. \end \end$$

(11)

The unconstrained OCP reads as in Eq. (9) except for the penalty function Eq. (11) instead of Eq. (8). As with penalty methods, a series of such are solved, but multiplier approximating function \(\lambda \) and penalty parameter \(\rho \) are updated according to specific rules [7].

Note that choosing \(\lambda = 0\) in Eq. (11) recovers Eq. (8). Thus, the penalty method, also known as Moreau–Yosida regularization in current context, can be regarded as a safeguarded augmented Lagrangian method with multiplier updates fixed to zero, see [7, Rem. 3.10]. A more elaborate multiplier update, i.e., approximating the Lagrange multiplier, might provide faster convergence [9] and improve numerical results.