Optimal control in pharmacokinetic drug administration

1.   Introduction

Optimal control theory has been applied to pharmacokinetic problems since the early 1980s: See the landmark papers ^[1] from 1983 and ^[2] from 1986. The first of the two mentioned papers studies in particular the optimal injection of a drug with respect to a certain optimization criterion. The second paper describes a global optimization method for finding the optimum of functions involving several real variables. As an application a method for the solution of certain optimal control problems associated to compartmental modelling was introduced.

In recent years interest in optimal control theory to solve problems of drug administration in terms of efficacy, toxicity, and overall costs has grown rapidly. See, e.g., the recent article ^[3] where the Control Theory for Therapy Design (CT4TD) framework is introduced, which employs optimal control theory on patient-specific pharmacokinetics and pharmacodynamics models, to deliver optimized therapeutic strategies. In ^[4] an optimal control problem for cancer chemotherapy is considered. The results are illustrated and discussed in the framework of a mathematical model for anti-angiogenic therapy. The authors of ^[5] consider optimization of drug release characteristics or dosing schedules for anticancer agents. They perform an optimal control analysis on a previously developed computational model for the testosterone effects of triptorelin in prostate cancer patients with the goal of finding optimal drug-release characteristics. In particular, numerical control optimization is used to find better therapeutic approaches in order to improve the final outcome of the patients. In the recent article ^[6] the authors investigate an optimal dosing strategy and develop an optimal dosing algorithm (OptiDose) that computes the optimal individualized dosing regimen for pharmacokinetic-pharmacodynamic models by using numerical methods. In ^[7] the Model Predictive Control (MPC) technology is discussed in the context of optimal drug administration. This method takes into account the constraints of the system and uses continuous plasma measurements to readjust the dosing schedule in an optimal and stable way at any time. Another recent approach is described in ^[8] where robust optimal control models are developed for optimal drug dosage design under model uncertainties.

With this article we address on the one hand experts in pharmacology who are interested in the mathematical tools used to solve problems of optimal dosing schedules. For these readers, we provide a brief and easily accessible introduction to the mathematical techniques. On the other hand, we would like to introduce mathematicians to interesting applications in the field of pharmacokinetic research. For them, the pharmacological context is briefly explained. Depending on the background, the reader can skip individual sections.

While past research has often focused on theoretical foundation, like existence and uniqueness of optimal control solutions, or numerical methods, the aim of this article is to provide a concrete way to actually compute explicit solutions for a class of optimal control problems occurring in drug administration. We keep the model as simple as possible so that it is easily accessible, while still keeping the framework sufficiently general to be applicable in real situations. To make the article self-contained we give a brief introduction to the theory of optimal control in Section 2. Section 3 discusses the application of optimal control in a pharmacokinetic two-box model. The optimal control problem for the central compartment is solved in Section 3.1. Finally, Section 3.2 deals with the corresponding problem for the peripheral compartment, which is subject to a natural non-linear constraint: Here we propose two approximate solutions which are easy to compute and to implement in the clinical context, and the exact solution by applying Pontryagin's Minimum Principle. At this point, we would like to refer the interested reader to the general literature on the theory of optimal control and the Pontryagin Minimum Principle, e.g., ^[9,10,11].

2.   The optimal control problem

We start with a general description of the problem of optimal control. An optimal control problem has two components:

1) The first basic component is a control system. Typically, such a control system is described by a system of ordinary differential equations:

Here $x$ is the state of the system that takes values in $\mathbb{R}^n$. Consider, for example, a compartmental model where the components of the vector $x$ indicate the concentration of an active pharmaceutical ingredient in the individual compartments. As usual, $t$ denotes the time. At time $t_0$ the initial state of the system is $x_0\in\mathbb R^n$. Additionally we have a control-input $u(t)$ depending on $t$ (also called control function or simply control) which we can influence, i.e., control. This is, e.g., the rate at which a drug is delivered to a compartment from the outside. Another example would be that $x$ describes the number of certain populations of an ecosystem and $u$ is an external influence, such as the input of pollutants, supply of food or the allowed shooting rate for hunting.

2) The second basic component is the cost functional. It assigns a "price" to each possible trajectory $x(t)$ of the system. Let us think of a control $u(t)$ given concretely and the system (1) solved for the given initial data. Let the corresponding "price" then be given by an integral

where $L$ is a given function and $T > 0$ is a (fixed in our case) time duration. That is, given the control $u(t)$, the resulting price $J(u)$ can be determined. Here, "price" usually does not mean a monetary amount, but simply a value to be optimized. For example, one would like that in a given compartment (organ) the solution (drug concentration) there deviates as little as possible from a given optimal target drug level during the time interval $(0, T)$: the average deviation (the "price") would thus have to be minimized. Or, in the case of the ecosystem, one would like to increase the population of an endangered species as much as possible: The sum of the animals of the population would be to maximize here.

This makes it already clear what the optimal control problem is all about:

Definition. The optimal control problem given by (1) and (2) is to find the control function $u(t)$ for which the cost functional $J(u)$ is minimum or maximum.

3.   Optimal control in a pharmacokinetic context

We illustrate the optimal control problem and its solution in an explicit example. For the general theoretical background we refer to the literature mentioned in the introduction.

Consider a compartment system with a central and a peripheral component: blood and an organ (see, e.g., [12,Chapter 6] for the use of compartment systems in pharmacokinetics). A drug is added to the blood at a certain rate $u(t)$: This will play the role of the control function here. In the blood and in the organ (for example the liver) the drug is degraded, while the two compartments exchange the substance at given rates. The dynamics of the process and the rates are described in Figure 1.

Let the concentrations of the drug in the blood and in the organ be denoted by $c_{\rm{B}}(t)$ and $c_{\rm{O}}(t)$, respectively. The corresponding system of differential equations for the vector

is given by

with the matrix

and the vector

At time $t = 0$ we suppose that both initial concentrations are zero: $c_{\rm{B}}(0) = c_{\rm{O}}(0) = 0$. In order to solve the system (3) we apply the Laplace transformation to it: Let

denote the Laplace transforms of the functions $c_{\rm{B}}, c_{\rm{O}}$ and $u$. Then the transformed system is

with

The solution of (3) is thus given by

where $\mathbb I$ denotes the identity matrix. Hence we have

Using partial fraction decomposition we write the solutions as

where $\lambda_1, \lambda_2$ are the eigenvalues^* of $A$, and

^*$\lambda_1+\lambda_2 = \operatorname{tr}A, \lambda_1\lambda_2 = \det A$. We only consider the case $\lambda_1\neq\lambda_2$ here. Observe that $\lambda_1, \lambda_2 < 0$.

In the form (6) the inverse Laplace transform is easy to carry out and by applying the appropriate rules one finds as solution of (3)

For a nice concrete example, we recommend the numerical values

in suitable units, which gives the integer eigenvalues $\lambda_1 = -10, \lambda_2 = -2$. We use these values to draw graphs of the solutions in the figures below.

3.1.   Optimal control in the central compartment

Suppose, the manufacturer of the drug that we use has defined an optimal concentration in the central compartment (blood) with maximum therapeutic effect, let us say a certain value $c_{\rm{B, opt}}$. It is therefore necessary to administer the drug in such a way so that during the treatment period $(0, T)$ the concentration $c_{\rm{B}}$ remains as closely as possible to this value $c_{\rm{B, opt}}$. This requirement can be formulated as follows:

Optimal control problem for the central compartment: The integral

is a measure for the deviation of the concentration $c_{\rm{B}}$ from the optimal value $c_{\rm{B, opt}}$ and should therefore be minimal!

The control function $u(t)$, namely the rate at which the drug is added to the blood, is present in $J(u)$ in the solution $c_{\rm{B}}$ which we found above in (7). If we now insert this solution in (10), we obtain

This integral has to be minimized by the optimal choice of the control function $u(t)$. For this purpose we use the calculus of variations: The basic idea of Leonhard Euler is to perturb the control function $u(t)$ a little bit. That is, we put in $J$ instead of $u$ the perturbed function $u+\varepsilon \varphi$. Here $\varepsilon$ is a real number, and $\varphi$ is a perturbation function that we can choose arbitrarily. Now we consider the function

If we assume that $u$ is the sought optimal control function that minimizes $J$, then $f$ has a local minimum at $0$, that is, $f'(0) = 0$ holds regardless of the choice of $\varphi$. If we calculate $f'(0)$ concretely, we get

The integrand of the outer integral is a product of two bracket expressions. In the first bracket, we are free to choose $\varphi$ as we wish. In particular, it can be shown that the first bracket can be made to represent any triangular function $d_{a, b}$ by appropriate choice of $\varphi$ (see Figure 2). Indeed this can even be explicitly achieved by setting the red bracket equal to $d_{a, b}$ and applying the Laplace transform.

In this way we can replace the first bracket expression in (12) by any triangular function $d_{a, b}$, and the integral must always yield the value $0$. The fundamental lemma of the calculus of variations now states that this is only possible if the function in the second bracket is identically equal to $0$ on the interval $(0, T)$: Namely, assume by way of contradiction that the second bracket is different from 0 in a point $t\in(0, T)$. Then, by continuity, it is strictly positive or strictly negative on an interval around $t$. But by choosing the support $(a-b, a+b)$ of $d_{ab}$ inside said interval, we get a value $f'(0)\neq 0$, contrary to our assumption. This argument shows that the optimal control function $u$ must satisfy the condition

for all $t\in(0, T)$. Observe that this condition actually means $J(u) = 0$. If we subject this condition to the Laplace transform, we get by (5)

We solve this equation for $U(s)$ and get by partial fraction decomposition

So now we only have to reconstruct the control function $u$ from its Laplace transform $U$. The original functions of the two fractions are easy to find:

However, observe that the constant function 1 does not correspond to an original function of the Laplace transform. But we can specify a sequence of functions whose Laplace transforms approximate arbitrarily well the function 1. Namely, we consider the sequence of functions

where $H$ denotes the Heaviside function: see Figure 3.

With this choice we indeed have

Here, we have used de L'Hôpital's rule for the calculation of the limit. For the formal limit of the sequence $\delta_n$ we write

$\delta$ is not a function but the Dirac delta distribution. We briefly comment on its interpretation for readers who are less familiar with this concept: Note that for any $n$ we always have

When an elastic ball hits the ground and bounces back, a large force acts on the ball in a short period of time. This force $F$ can therefore be expressed (in suitable units) by a function $F(t) = \delta_n(t)$. The integral $\int_{-\infty}^\infty F(t)dt = 1$ is then just the (constant) change in momentum of the ball, which it experiences when bouncing back. The harder the ball, the shorter the period of contact with the ground. The distribution $\delta$ is therefore an idealization of such an impulse.

But let us go back to our optimal control problem. Putting together (14)–(16) we are now ready to state the following.

Theorem 1. The control function

solves the optimal control problem for the central compartment, and $J(u) = 0$.The limits are

and $u$ is monotone decreasing. In particular, $u$ satisfies automatically the natural constraint $u(t)\ge 0$.

Clinical interpretation of the result. Optimal treatment means that the patient is given an injection at the beginning $t = 0$ as a bolus with the amount that would, if it were not degraded, raise the concentration of the drug in the blood to the optimal concentration $c_{\rm{B, opt}}$: This corresponds to the term $c_{\rm{B, opt}}\delta$ (see, e.g., [12,Chapter 7] for the pharmacokinetics of an intravenous bolus injection). At the same time, the infusion (drip) is set up so that the delivery occurs according to $\frac{c_{\rm{B, opt}}}{k_{\rm O} + k_{\rm{BO}}} \Bigl(\det A+ e^{-(k_{\rm O} + k_{\rm{BO}})t} k_{\rm{BO}} k_{\rm{OB}} \Bigr)$. Such dosages are possible via programmable infusion pumps (see, e.g., [12,Chapter 11]).

The left hand side of Figure 4 shows the solution $c_{\rm{B}}$, if $\delta$ is approximated by the sequence $\delta_n$: Also in reality, the administration of a syringe requires a certain time, i.e., a time interval $(0, 1/n)$.

Somewhat disturbing is the overshooting of the optimal concentration. The reason is that the infusion starts at the same time as the bolus, so the two doses add up in the time interval $(0, 1/n)$ (see Figure 4 on the right). In an extreme case and for a drug with a very narrow therapeutic window, the overshoot might even reach a toxic level. Although the excess becomes smaller when the time period for the bolus is shortened, the overshoot remains. This can be avoided if the infusion is switched on only after the administration of the syringe, i.e., with the control function $u(t) = \frac{c_{\rm{B, opt}}}{k_{\rm O}+k_{\rm{BO}}}(\det A+ e^{-(k_{\rm O}+k_{\rm{BO}})(t-1/n)} k_{\rm{BO}}k_{\rm{OB}})H(t-\frac1n)$. Figure 5 shows on the left the behaviour of the solution $c_{\rm{B}}$ in the case of such a delayed infusion, and the corresponding control function on the right.

3.2.   Optimal control in the peripheral compartment

The optimal control problem for the peripheral compartment (the organ, e.g., the liver) is significantly more delicate for the following reason: The control function $u$ which we found for the central compartment respects automatically the natural side condition $u\ge 0$. This is not the case for the peripheral compartment as we shall see below. We therefore have to implement the non-linear constraint $u\ge0$ explicitly in the formulation of the problem:

Optimal control problem for the peripheral compartment: Minimize the deviation of the concentration $c_{\rm{O}}$ from the optimal therapeutic concentration $c_{\rm{O, opt}}$, i.e., the cost functional

under the constraint $u\ge0$.

With the solution for $c_{\rm{O}}$ in (8) we obtain the explicit expression

Neglecting the constraint, the optimality condition which corresponds to (13) turns out to be

According to (5) the Laplace transformed version of this equation is

which gives

Applying the inverse Laplace transformation leads to

This solution, however, violates the non-linear, pointwise constraint $u(t)\ge 0$ for any reasonable approximation of the delta distribution. E.g., we have $\delta_n'(t) = n\bigl(\delta(t)-\delta(t-1/n)\bigr)$. Indeed, pushing the level of the concentration in the peripheral compartment to the desired level $c_{\rm{O, opt}}$ as fast as possible would require a massive bolus at time $t = 0$. But then one would have to reduce (by negative values of $u$) the concentration in the central compartment in order to avoid that the flow from the central to the peripheral compartment raises the concentration there far beyond the level $c_{\rm{O, opt}}$.

Before we address the problem of this constraint, it is instructive to consider the dual problem: Observe that the control function $u$ determines the solution $c_{\rm{O}}$ in (8). But vice versa, a given function $c_{\rm{O}}$ determines the corresponding control function $u$: Indeed, by (5) we have

The inverse Laplace transformation yields

where we have used that $c_{\rm{O}}(0) = 0$. We can therefore reformulate the optimal control problem as follows:

Dual formulation of the optimal control problem for the peripheral compartment: The integral

is to be minimized under the side conditions

In this formulation, the control function $u$ has disappeared. It is no longer a problem of optimal control, but simply a problem in the calculus of variations with inequality constraints (such problems can be treated by the method of slack variables: see, e.g., ^[13] or ^[14]). We use it in the following sections because (20) indicates how the solution $c_{\rm O}$ behaves where the differential inequality is saturated.

3.2.1.   A dosage scheme without overshoot

In view of the dual formulation of the optimal control problem, a reasonable approach is to choose $c_{\rm{O}}(t)$ as a $C^1$ function with $c_{\rm{O}}(t)\equiv c_{\rm{O, opt}}$ on an interval $(T_0, T]$ such that the integrand in (19) vanishes on this interval, and $c_{\rm{O}}(t)\det A -c'_{\rm{O}}(t)\operatorname{tr} A+c''_{\rm{O}}(t)\equiv0$ on $[0, T_0]$, with $c_{\rm{O}}(0) = 0$ (compare (20)). Hence, $c_{\rm{O}}(t)$ has the form

In order to obtain a global $C^1$ function we have to choose $\varkappa$ and $T_0$ in (21) such that $c_{\rm{O}}(T_0) = c_{\rm{O, opt}}$ and $c'_{\rm{O}}(T_0) = 0$. A short calculation gives

So, we obtain for $c_{\rm{O}}$ and the corresponding control function $u$, given by (18),

(see Figure 6).

Clinical interpretation. The dosing scheme in (23) means that at time $t = 0$ a bolus with a dose $-\frac{c_{\rm{O, opt}}}{k_{\rm{OB}}} \lambda_1 \bigl(\frac{\lambda_1}{\lambda_2}\bigr)^{\lambda_2/(\lambda_1-\lambda_2)}$ is given^†, and at time $t = T_0$ a constant drip $\frac{c_{\rm{O, opt}}\det A}{k_{\rm{OB}}}$ per time unit starts.

^†recall, that $\lambda_1, \lambda_2 < 0$

The nice feature of this dosing scheme is, that its parameters are explicitly given in terms of the constants governing the pharmacokinetic process, and that the optimal concentration $c_{\rm{O, opt}}$ is never exceeded. Moreover, it requires only one syringe at time $t = 0$. However, this is not the solution of the optimal control problem, since $J$ is not minimal with this choice of $u$.

3.2.2.   Restricting the class of control functions

If we allow a mild overshoot of the concentration in the organ, we can slightly improve the dosing schedule from Section 3.2.1. To do this we restrict the set of the control functions to a finite dimensional space: We minimize the original cost functional (17) for the class $\mathcal C$ of piecewise constant functions of the form

We accept here that the minimum of the cost function in this smaller class $\mathcal C$ is larger than if we allow all admissible control functions. The advantage is, on the other hand, that we have to determine the minimum only in a finite dimensional set.

In order to find the optimal values $\bar \varkappa, \bar b$ and $\bar T_0$ we proceed as follows: According to (8) we obtain for a control function $u$ in the class $\mathcal C$

Since the therapy is likely to take a long time, we set $T = \infty$ in (17). Then, for the integral in (17) to be finite, $\lim_{t\to\infty} c_{\rm{O}} (t) = c_{\rm{O, opt}}$ must hold. This is the case if and only if $b = \bar b = \frac{c_{\rm{O, opt}}\det A}{k_{\rm{OB}}}$. Then, $J(u)$ in (17) becomes a function in the variables $\varkappa$ and $T_0$:

This function has on the set $\{(\varkappa, T_0) : \varkappa > 0, T_0 > 0\}$ a global minimum where its gradient

vanishes. The numerical solution of $\nabla J = 0$ then yields the optimal values for $\varkappa$ and $T_0$. Hence, we have:

Theorem 2. The optimal dosing schedule among the control functions in the class $\mathcal C$ for the peripheral compartment is$u(t) = \bar\varkappa\delta(t) + \frac{c_{\mathit{\rm{O, opt}}}\det A}{k_{\rm{OB}}} H(t-\bar T_0)$ where $\bar\varkappa$ and $\bar T_0$ are thepositive solutions of the system

The corresponding graphs are shown in Figure 7.

Clinical interpretation. Compared to the values in (22), the initial bolus is slightly larger, and the waiting time for the drip to start a bit longer. For the numerical values in (9), the value of $J$ for this control function is about 9% lower than for the control function in (23).

3.2.3.   Solution with Pontryagin's Minimum Principle

We now apply Pontryagin's Minimum Principle in order to solve the optimal control problem in the peripheral compartment exactly. To do this, we consider the Hamiltonian of the problem: In general, if we want to minimize $J = \int_0^TL(x(t), u(t))dt$ where $x$ solves the system $\dot x = F(x, u), \ x(0) = x_0$, with control function $u$, the Hamiltonian is given by $H(x, u, p) = L(x, u)+p^\top F(x, u)$. Here, $p$ is a Lagrange multiplier, and now only control functions $u$ in a given class $U$ are allowed. Suppose $\bar x$ is the solution for the optimal control function $\bar u$ and the corresponding Lagrange multiplier $\bar p$. Then Pontryagin's Minimum Principle states that $\bar u$ minimizes $H$ for all $t\in(0, T)$, i.e.,

holds for all admissible control functions $u\in U$. In addition the state and co-state equations

must hold.

In our case we have

where $a_{ij}$ are the coefficients of the matrix $A$. The class $U$ of admissible control functions consists of all functions $u$ for which $u(t)\ge 0$ for all $t\in[0, T]$. So, we have to solve the following state and co-state equations:

Observe that only the term $p_1u$ in the Hamiltonian depends on $u$. Therefore the condition (25) for the optimal control function $\bar u$ reduces to the following:

This means, that on an interval where $\bar p_1 > 0$ the control function $\bar u$ must vanish identically, and $\bar u(t) > 0$ is only possible where $\bar p_1(t)\le 0$. Moreover, $\bar p_1 < 0$ cannot hold on a set of positive measure, since there $u = \infty$ would lead to $J = \infty$. Moreover, observe that on an interval where $p_1$ vanishes, also $p_2$ vanishes by (28) and therefore, by (29), $c_{\rm O}$ must be identical to the optimal value $c_{\rm{O, opt}}$ on said interval.

For simplicity, we assume again that $T = \infty$. The idea to construct the solution to (26)–(30) is the following: Suppose $p_1$ vanishes on a maximum interval $[0, \infty)$, then we construct the solution backward in time, and finally we shift everything in time to the right so that the zero of $c_{\rm O}$ is at $t = 0$. In concrete terms, we proceed as follows: Let $X(t) = (f(t), g(t))^\top$ be a solution of the homogeneous system $\dot X = AX$ with $g(0) = c_{\rm{O, opt}}$. Then, $g$ is of the form

where $a$ is a free parameter which will be determined below. Observe that $g_a$ satisfies $g_a(t)\det A-g'_a(t)\operatorname{tr}A+g''_a(t) = 0$ (compare (20)). Both eigenvalues of $A$ are strictly negative, and we assume that $\lambda_1 < \lambda_2 < 0$. If the parameter $a$ satisfies

then we have

Then one solves the systems (28) and (29) with $p_1(0) = p_2(0) = 0$ with $g_a$ in place of $c_{\rm O}$ in (29). The resulting solution $p_1, p_2$ depends on $a$ and can be written down explicitly, but the expression is unpleasantly long. It turns out, however, that $p_1$ is strictly positive on a maximum interval $(T_1(a), 0)$ with $T_1(a) < 0$ depending on $a$. Now, choose the parameter $a = \bar a$ such that $T_1(\bar a) = T_0(\bar a)$: For this value $a = \bar a$, the functions $g_{\bar a}$ and $p_1$ are both strictly positive on $(T_0, 0)$ and vanish in $T_0$. How do we find this value $\bar a$? The function $a\mapsto p_1(T_0(a))$ has a negative minimum in $a = \frac{\lambda_2}{\lambda_2-\lambda_1} < 0$ and is monotone increasing with $\lim_{a\nearrow 0}p_1(T_0(a)) = \infty$. Therefore the equation $p_1(T_0(a)) = 0$ has a unique solution $\bar a\in(\frac{\lambda_2}{\lambda_2-\lambda_1}, 0)$ which can be determined numerically.

Now we are almost done. We only have to shift the constructed functions to the right, so that the zero $T_0$ ends up at $t = 0$:

Theorem 3. With the notation above, the optimal control solution for the peripheral compartment is

and the control function $u$ is given by

Indeed, this solution satisfies the conditions of the Pontryagin Minumum Principle by construction. Observe that $u$ is given by (18): It is important to note that $c'_{\rm{O}}$ is discontinuous in $T_0$. Therefore, in (18) we have $c''_{\rm{O}}(T_0) = - g'_{\bar a}(0)\delta(t-T_0)$. Notice also that $c_{\rm B}$ follows easily from (27).

Clinical interpretation of the result. The optimal treatment is as follows:

${\rm{i)}}$ A bolus with the dose $\frac{c_{\rm{O}}'(0)}{k_{\rm{OB}}}$ at time $t = 0$

${\rm{ii)}}$ No infusion in the time interval $(0, T_0)$

${\rm{iii)}}$ A bolus with the dose $-\frac{g_{\bar a}'(0)}{k_{\rm{OB}}}$ at time $t = T_0$

${\rm{iv)}}$ A constant drip with $\frac{c_{\rm{O, opt}}\det A}{k_{\rm{OB}}}$ per time unit for $t > T_0$.

Figure 8 shows the solution for the numerical values in (9). In addition, Figure 9 shows the corresponding solution $c_{\rm B}$ and the co-state solution $p_1$. See also [12,Chapter 12] for multiple intravenous bolus injections.

We close the discussion with two final remarks:

● A posteriori, the optimal control solution which we found by the Pontryagin Minimum Principle is optimal not only for $T = \infty$, but for all $T\ge T_0$.

● For the numerical values in (9) the cost functional $J$ is only about 0.5% smaller for the Pontryagin solution compared to the solution which we found in Section 3.2.2.

4.   General conclusions

The optimal control problem for the central compartment is easy to solve explicitly (Section 3.1). In particular, the natural constraint $u\ge 0$ is automatically satisfied. In contrast, this constraint must be explicitly considered for the solution in the peripheral compartment. If one is only interested in an approximate solution, different variants are available: An approximate solution without overshoot can easily be specified explicitly in terms of the parameters of the system (Section 3.2.1). Numerically, it is also quite easy to calculate the parameters of an approximate solution if one restricts the control functions to a finite-dimensional space (Section 3.2.2). With somewhat greater numerical effort, the parameters of the exact solution of the optimal control problem can be determined using the Pontryagin Minimum Principle (Section 3.2.3). In each case, the dosage schedule can be easily derived from the respective result. The two approximate solutions get by with only one syringe, while the Pontryagin solution requires two syringes.

Acknowledgements

I would like to thank the referees for their careful reading and the valuable suggestions which significantly helped to improve the content of the article. I also thank Chris Busenhart and his father for the reference to the work of Hans P. Geering.

Conflict of interest

The author declares there is no conflict of interest.