An Analytic Method to Determine the Optimal Time for the Induction Phase of Anesthesia

We obtain an analytical solution for the time-optimal control problem in the induction phase of anesthesia. Our solution is shown to align numerically with the results obtained from the conventional shooting method. The induction phase of anesthesia relies on a pharmacokinetic/pharmacodynamic (PK/PD) model proposed by Bailey and Haddad in 2005 to regulate the infusion of propofol. In order to evaluate our approach and compare it with existing results in the literature, we examine a minimum-time problem for anesthetizing a patient. By applying the Pontryagin minimum principle, we introduce the shooting method as a means to solve the problem at hand. Additionally, we conducted numerical simulations using the MATLAB computing environment. We solve the time-optimal control problem using our newly proposed analytical method and discover that the optimal continuous infusion rate of the anesthetic and the minimum required time for transition from the awake state to an anesthetized state exhibit similarity between the two methods. However, the advantage of our new analytic method lies in its independence from unknown initial conditions for the adjoint variables.


Introduction
Based on Guedel's classification, the first stage of anesthesia is the induction phase, which begins with the initial administration of anesthesia and ends with loss of consciousness [1].Millions of people safely receive several types of anesthesia while undergoing medical procedures: local anesthesia, regional anesthesia, general anesthesia, and sedation [2].However, there may be some potential complications of anesthesia including anesthetic awareness, collapsed lung, malignant hyperthermia, nerve damage, and postoperative delirium.Certain factors make it riskier to receive anesthesia, including advanced age, diabetes, kidney disease, heart disease, high blood pressure, and smoking [3].To avoid the risk, administering anesthesia should be carried out on a scientific basis, based on modern pharmacotherapy, which relies on both pharmacokinetic (PK) and pharmacodynamic (PD) information [4].Pharmacokinetics is used to describe the absorption and distribution of anesthesia in body fluids, resulting from the administration of a certain anesthesia dose.Pharmacodynamics is the study of the effect resulting from anesthesia [5].Multiple mathematical models were already presented to predict the dynamics of the pharmacokinetics/pharmacodynamics (PK/PD) models [6][7][8][9].Some of these models were implemented following different methods [2,10,11].
The parameters of PK/PD models were fitted by Schnider et al. in [12].In [6], the authors study pharmacokinetic models for propofol, comparing Schnider et al. and Marsh et al. models [13].The authors of [6] conclude that Schnider's model should always be used in effect-site targeting mode, in which larger initial doses are administered but smaller than those obtained from Marsh's model.However, users of the Schnider model should be aware that in the morbidly obese, the lean body mass (LBM) equation can generate paradoxical values, resulting in excessive increases in maintenance infusion rates [12].In [14], a new strategy is presented to develop a robust control of anesthesia for the maintenance phase, taking into account the saturation of the actuator.The authors of [15] address the problem of optimal control of the induction phase.For other related works, see [8,16] and references therein.
Here, we consider the problem proposed in [15], to transfer a patient from a state consciousness to unconsciousness.We apply the shooting method [17] using the Pontryagin minimum principle [18], correcting some inconsistencies found in [15] related with the stop criteria of the algorithm and the numerical computation of the equilibrium point.Secondly, we provide a new different analytical method to the time-optimal control problem for the induction phase of anesthesia.While the shooting method, popularized by Zabi et al. [15], is widely employed for solving such control problems and determining the minimum time, its reliance on Newton's method makes it sensitive to initial conditions.The shooting method's convergence is heavily dependent on the careful selection of initial values, particularly for the adjoint vectors.To overcome this limitation, we propose an alternative approach, which eliminates the need for initial value selection and convergence analysis.Our method offers a solution to the time-optimal control problem for the induction phase of anesthesia, free from the drawbacks associated with the shooting method.Furthermore, we propose that our method can be extended to other PK/PD models to determine optimal timings for drug administration.To compare the methods, we perform numerical simulations to compute the minimum time to anesthetize a man of 53 years, 77 kg, and 177 cm, as considered in [15].We find the optimal continuous infusion rate of the anesthetic and the minimum time that needs to be chosen for treatment, showing that both the shooting method of [15] and the one proposed here coincide.
This paper is organized as follows.In Section 2, we recall the pharmacokinetic and pharmacodynamic model of Bailey and Haddad [19], the Schnider model [12], the bispectral index (BIS), and the equilibrium point [14].Then, in Section 3, a time-optimal control problem for the induction phase of anesthesia is posed and solved both by the shooting and analytical methods.Finally, in Section 4, we compute the parameters of the model using the Schnider model [12], and we illustrate the results of the time-optimal control problem through numerical simulations.We conclude that the optimal continuous infusion rate for anesthesia and the minimum time that should be chosen for this treatment can be found by both shooting and analytical methods.The advantage of the new method proposed here is that it does not depend on the concrete initial conditions, while the shooting method is very sensitive to the choice of the initial conditions of the state and adjoint variables.We end with Section 5 of conclusions, pointing also some directions for future research.

The PK/PD Model
The pharmacokinetic/pharmacodynamic (PK/PD) model consists of four compartments: intravascular blood (x 1 (t)), muscle (x 2 (t)), fat (x 3 (t)), and effect site (x 4 (t)).The effect site compartment (brain) is introduced to account for the finite equilibration time between central compartment and central nervous system concentrations [19].This model is used to describe the circulation of drugs in a patient's body, being expressed by a four-dimensional dynamical system as follows: x 1 (t) − a e 0 x 4 (t). (1) The state variables for system (1) are subject to the following initial conditions: where x 1 (t), x 2 (t), x 3 (t), and x 4 (t) represent, respectively, the masses of the propofol in the compartments of blood, muscle, fat, and effect site at time t.The control u(t) is the continuous infusion rate of the anesthetic.The parameters a 1 0 and a e 0 represent, respectively, the rate of clearance from the central compartment and the effect site.The parameters a 1 2 , a 1 3 , a 2 1 , a 3 1 , and a e 0 /v 1 are the transfer rates of the drug between compartments.A schematic diagram of the dynamical control system (1) is given in Figure 1.

Schnider's Model
Following Schnider et al. [12], the lean body mass (LBM) is calculated using the James formula, which performs satisfactorily in normal and moderately obese patients, but not so well for severely obese cases [20].The James formula calculates LBM as follows: for Female, LBM = 1.07 The parameters of the PK/PD model ( 1) are then estimated according to Table 1.

The Bispectral Index (BIS)
The BIS is the depth of anesthesia indicator, which is a signal derived from the EEG analysis and directly related to the effect site concentration of x 4 (t).It quantifies the level of consciousness of a patient from 0 (no cerebral activity) to 100 (fully awake patient), and can be described empirically by a decreasing sigmoid function [19]: where BIS 0 is the BIS value of an awake patient typically set to 100, EC 50 corresponds to the drug concentration associated with 50% of the maximum effect, and γ is a parameter modeling the degree of nonlinearity.According to [21], typical values for these parameters are EC 50 = 3.4 mg/L and γ = 3.

The Equilibrium Point
Following [14], the equilibrium point is obtained by equating the right-hand side of (1) to zero, with the condition It results that the equilibrium point x e = (x e 1 , x e 2 , x e 3 , x e 4 ) is given by and the value of the continuous infusion rate for this equilibrium is The fast state is defined by The control of the fast dynamics is crucial because the BIS is a direct function of the concentration at the effect site.

Time-Optimal Control Problem
Let We can write the dynamical system (1) in a matrix form as follows: where Here, the continuous infusion rate u(t) is to be chosen so as to transfer the system (1) from the initial state (wake state) to the fast final state (anesthetized state) in the shortest possible time.Mathematically, we have the following time-optimal control problem [15]: where t f is the first instant of time that the desired state is reached, and C and x eF are given by with

Pontryagin Minimum Principle
According to the Pontryagin minimum principle (PMP) [18], if ũ ∈ L 1 is optimal for Problem (13) and the final time t f is free, then there exists where the Hamiltonian H is defined by Moreover, the minimality condition holds almost everywhere on t ∈ [0, t f ].
Since the final time t f is free, according to the transversality condition of PMP, we obtain: Solving the minimality condition (18) on the interior of the set of admissible controls gives the necessary condition where ψ1 (t) is obtained from the adjoint system (16), that is, ψ′ (t) = −A T ψ(t), and the transversality condition (19).This is discussed in Sections 3.2 and 3.3.

Shooting Method
The shooting method is a numerical technique used to solve boundary value problems, specifically in the realm of differential equations and optimal control.It transforms the problem into an initial value problem by estimating the unknown boundary conditions.Through iterative adjustments to these estimates, the boundary conditions are gradually satisfied.In [17], the authors propose an algorithm that addresses numerical solutions for parameterized optimal control problems.This algorithm incorporates multiple shooting and recursive quadratic programming, introducing a condensing algorithm for linearly constrained quadratic subproblems and high-rank update procedures.The algorithm's implementation leads to significant improvements in convergence behavior, computing time, and storage requirements.For more on numerical approaches to solve optimal control problems, we refer the reader to [22] and references therein.

Analytical Method
We now propose a different method to choose the optimal control.If the pair (A, B) satisfies the Kalman condition and all eigenvalues of matrix A ∈ n × n are real, then any extremal control has at most n − 1 commutations on R + (at most n − 1 switching times).We consider the following eight possible strategies: Strategy 1 (zero switching times): Strategy 2 (zero switching times): Strategy 3 (one switching time): where t c is a switching time.
Strategy 4 (one switching time): Strategy 5 (two switching times): where t c1 and t c2 represent two switching times.Strategy 6 (two switching times): Strategy 7 (three switching times): where t c1 , t c2 , and t c3 represent three switching times.
Strategy 8 (three switching times): Let x(t) be the trajectory associated with the control u(t), given by the relation where exp(A) is the exponential matrix of A.
To calculate the switching times t c , t c1 , t c2 and the final time t f , we have to solve the following nonlinear equation: xeF (t f ) = (x e1 , x e4 ). (37) We also solve (37) using the Newton method [23].

Numerical Example
In this section, we use the shooting and analytical methods to calculate the minimum time t f to anesthetize a man of 53 years, 77 kg, and 177 cm.
The equilibrium point and the flow rate corresponding to a BIS of 50 are: x e = (14.518mg, 64.2371 mg, 813.008 mg, 3.4 mg), u e = 6.0907 mg/min.( 38) Following the Schnider model, the matrix A of the dynamic system ( 11) is given by: We are interested to solve the following minimum-time control problem: (40)

Numerical Resolution by the Shooting Method
Let z(t) = (x(t), ψ(t)).We consider the following Cauchy problem: where The shooting function S is given by where All computations were performed with the MATLAB numeric computing environment, version R2020b, using the medium-order method and the function ode45 (Runge-Kutta method) in order to solve the nonstiff differential system (22).We have used the variable order method and the function ode113 (Adams-Bashforth-Moulton method) in order to solve the nonstiff differential system (25), and the function fsolve in order to solve equation S(z 0 ) = 0. Thus, we obtain that the minimum time is equal to with ψ T (0) = (−0.0076,0.0031, −0.0393, −0.0374). (46)

Numerical Resolution by the Analytical Method
The pair (A, B) satisfies the Kalman condition, and the matrix A has four real eigenvalues.Then, the extremal control u(t) has at most three commutations on R + .Therefore, let us test the eight strategies provided in Section 3.3.
Note that the anesthesiologist begins with a bolus injection to transfer the patient state from the consciousness state x(0) to the unconsciousness state x eF = (14.518,3.4), that is, u(0) = U max = 106.0907mg/min.(47) Thus, Strategies 2, 4, 6, and 8 are not feasible here.Therefore, in the sequel, we investigate Strategies 1, 3, 5, and 7 only.
Strategy 1: Let u(t) = 106.0907mg/min for all t ∈ [0, t f ].The trajectory x(t), associated with this control u(t), is given by the following relation: System (37) takes the form and has no solutions.Thus, Strategy 1 is not feasible.
Strategy 3: Let u(t), t ∈ [0, t f ], be the control defined by The trajectory x(t) associated with this control u(t) is given by To calculate the switching time t c and the final time t f , we have to solve the nonlinear system (52) with the new condition Similarly to Section 4.1, all numerical computations were performed with MATLAB R2020b using the command solve to solve Equation (52).The obtained minimum time is equal to with the switching time t c = 0.5467 min.
Strategy 5: Let u(t), t ∈ [0, t f ], be the control defined by the relation where t c1 and t c2 are the two switching times.The trajectory x(t) associated with control (59) is given by To compute the two switching times t c1 and t c2 and the final time t f , we have to solve the nonlinear system (52) with 0 ≤ t c1 ≤ t c2 ≤ t f .(61) It turns out that System (52) subject to Condition (61) has no solution.Thus, Strategy 5 is also not feasible.
Strategy 7: Let u(t), t ∈ [0, t f ], be the control defined by the relation where t c1 , t c2 , and t c3 are the three switching times.The trajectory x(t) associated with Control (62) is given by To compute the three switching times t c1 , t c2 , and t c3 and the final time t f , we have to solve the nonlinear system (52) with It turns out that System (52) subject to Condition (64) has no solution.Thus, Strategy 7 is also not feasible.
In Figures 2 and 3, we present the solutions of the linear system of differential equations (40) under the optimal control u(t) illustrated in Figure 4, where the black curve corresponds to the one obtained by the shooting method, as explained in Section 3.2, while the blue curve corresponds to our analytical method, in the sense of Section 3.3.In addition, for both figures, we show the controlled BIS Index, the trajectory of fast states corresponding to the optimal continuous infusion rate of the anesthetic u(t), and the minimum time t f required to transition System (40) from the initial (wake) state x 0 = (0, 0, 0, 0) to the fast final (anesthetized) state x eF = (14.518,3.4) in the shortest possible time.The minimum time t f is equal to t f = 1.8397 min by the shooting method (black curve in Figure 2), and it is equal to t f = 1.8397 min by the analytical method (blue curve in Figure 3).
By using the shooting method, the black curve in Figure 4 shows that the optimal continuous infusion rate of the induction phase of anesthesia u(t) is equal to 106.0907 mg/min until the switching time t c = 0.5467 min.
Then, it is equal to 0 mg/min (stop-infusion) until the final time t f = 1.8397 min,  By using the analytical method, the blue curve in Figure 4 shows that the optimal continuous infusion rate of the induction phase of anesthesia u(t) is equal to 106.0907 mg/min until the switching time t c = 0.5467 min.
Then, it is equal to 0 mg/min (stop-infusion) until the final time t f = 1.8397 min.
We conclude that both methods work well and give similar results.However, in general, the shooting method does not always converge, depending on the initial conditions (46).To obtain such initial values is not an easy task since no theory is available to find them.For this reason, the proposed analytical method is logical, practical, and more suitable for real applications.

Conclusions
The approach proposed by the theory of optimal control is very effective.The shooting method was proposed by Zabi et al. [15], which is used to solve the time-optimal control problem and calculate the minimum time.However, this approach is based on Newton's method.The convergence of Newton's method depends on the initial conditions, being necessary to select an appropriate initial value so that the function is differentiable and the derivative does not vanish.This implies that the convergence of the shooting method is attached to the choice of the initial values.Therefore, the difficulty of the shooting method is to find the initial conditions of the adjoint vectors.Here, the aim was to propose a different approach, which we call "the analytical method", that allows to solve the timeoptimal control problem for the induction phase of anesthesia without such drawbacks.Our method is guided by the selection of the optimal strategy, without the need to choose initial values and study the convergence.We claim that our method can also be applied to other PK/PD models, in order to find the optimal time for the drug administration.
In the context of PK/PD modeling, the challenges associated with uncertainties in plant model parameters and controller gains for achieving robust stability and controller non-fragility are significant [24].These challenges arise from factors like inter-individual variability, measurement errors, and the dynamic nature of patient characteristics and drug response.Further investigation is needed to understand and develop effective strategies to mitigate the impact of these uncertainties in anesthesia-related PK/PD models.This research can lead to the development of robust and non-fragile control techniques that enhance the stability and performance of anesthesia delivery systems.By addressing these challenges, we can improve the precision and safety of drug administration during anesthesia procedures, ultimately benefiting patient outcomes and healthcare practices.In this direction, the recent results of [25] may be useful.Moreover, we plan to investigate PK/PD fractional-order models, which is a subject under strong current research [26].This is under investigation and will be addressed elsewhere.

Figure 1 .
Figure 1.Schematic diagram of the PK/PD model with the effect site compartment of Bailey and Haddad [19].

Figure 2 .
Figure 2. The state trajectory, controlled BIS index, and trajectory of the fast states corresponding to the optimal control u(t) of Figure 4, using the shooting method.

Figure 3 .Figure 4 .
Figure 3.The state trajectory, controlled BIS index, and trajectory of the fast states corresponding to the optimal control u(t) of Figure 4, using the analytical method.