Next Article in Journal
Zalcman Functional and Majorization Results for Certain Subfamilies of Holomorphic Functions
Next Article in Special Issue
Some Axioms and Identities of L-Moments from Logistic Distribution with Generalizations
Previous Article in Journal
Fractional p-Laplacian Coupled Systems with Multi-Point Boundary Conditions
Previous Article in Special Issue
Reliability Estimation of Inverse Weibull Distribution Based on Intuitionistic Fuzzy Lifetime Data
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

by
Mohamed A. Zaitri
1,2,
Cristiana J. Silva
1,3,* and
Delfim F. M. Torres
1,4
1
Center for Research and Development in Mathematics and Applications (CIDMA), Department of Mathematics, University of Aveiro, 3810-193 Aveiro, Portugal
2
Department of Mathematics, University of Djelfa, Djelfa 17000, Algeria
3
Department of Mathematics, ISCTE—Instituto Universitário de Lisboa, 1649-026 Lisbon, Portugal
4
Research Center in Exact Sciences (CICE), Faculty of Sciences and Technology (FCT), University of Cape Verde (Uni-CV), Praia 7943-010, Cape Verde
*
Author to whom correspondence should be addressed.
Axioms 2023, 12(9), 867; https://doi.org/10.3390/axioms12090867
Submission received: 10 June 2023 / Revised: 4 September 2023 / Accepted: 6 September 2023 / Published: 8 September 2023
(This article belongs to the Special Issue Mathematical Methods in the Applied Sciences)

Abstract

:
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.

1. 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.

2. 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 1 0 + a 1 2 + a 1 3 ) x 1 ( t ) + a 2 1 x 2 ( t ) + a 3 1 x 3 ( t ) + u ( t ) , x ˙ 2 ( t ) = a 1 2 x 1 ( t ) a 2 1 x 2 ( t ) , x ˙ 3 ( t ) = a 1 3 x 1 ( t ) a 3 1 x 3 ( t ) , x ˙ 4 ( t ) = a e 0 v 1 x 1 ( t ) a e 0 x 4 ( t ) .
The state variables for system (1) are subject to the following initial conditions:
x ( 0 ) = x 1 ( 0 ) , x 2 ( 0 ) , x 3 ( 0 ) , x 4 ( 0 ) = 0 , 0 , 0 , 0 ,
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.

2.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 Male , LBM = 1.1 × weight 128 × weight height 2 ,
for Female , LBM = 1.07 × weight 148 × weight height 2 .
The parameters of the PK/PD model (1) are then estimated according to Table 1.

2.2. 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]:
B I S ( x 4 ( t ) ) = B I S 0 1 x 4 ( t ) γ x 4 ( t ) γ + E C 50 γ ,
where B I S 0 is the B I S value of an awake patient typically set to 100, E C 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 E C 50 = 3.4 mg/L and γ = 3 .

2.3. The Equilibrium Point

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

3. Time-Optimal Control Problem

Let x ( t ) = ( x 1 ( t ) , x 2 ( t ) , x 3 ( t ) , x 4 ( t ) ) R 4 . We can write the dynamical system (1) in a matrix form as follows:
x ˙ ( t ) = A x ( t ) + B u ( t ) ,
where
A = ( a 1 0 + a 1 2 + a 1 3 ) a 2 1 a 3 1 0 a 1 2 a 2 1 0 0 a 1 3 0 a 3 1 0 a e 0 v 1 0 0 a e 0 and B = 1 0 0 0 .
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]:
min t f J = 0 t f d t , x ˙ ( t ) = A x ( t ) + B u ( t ) , x ( 0 ) = ( 0 , 0 , 0 , 0 ) , C x e F ( t f ) = x e F , 0 u ( t ) U m a x , t [ 0 , t f ] , t f is free ,
where t f is the first instant of time that the desired state is reached, and C and x e F are given by
C = 1 0 0 1 , x e F = ( x e 1 , x e 4 ) ,
with
x e F ( t f ) = ( x 1 ( t f ) , x 2 ( t f ) ) .

3.1. Pontryagin Minimum Principle

According to the Pontryagin minimum principle (PMP) [18], if u ˜ L 1 is optimal for Problem (13) and the final time t f is free, then there exists ψ ( t ) = ( ψ 1 ( t ) , , ψ 4 ( t ) ) , t [ 0 , t f ] , ψ A C ( [ 0 , t f ] ; R 4 ) , called the adjoint vector, such that
x ˙ = H ψ , ψ ˙ = H x ,
where the Hamiltonian H is defined by
H ( t , x , u , ψ ) = 1 + ψ T ( A x + B u ) .
Moreover, the minimality condition
H ( t , x ˜ ( t ) , u ˜ ( t ) , ψ ˜ ( t ) ) = min 0 u U m a x H ( t , x ˜ ( t ) , u , ψ ˜ ( t ) )
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:
H ( t f , x ( t f ) , u ( t f ) , ψ ( t f ) ) = 0 .
Solving the minimality condition (18) on the interior of the set of admissible controls gives the necessary condition
u ˜ ( t ) = 0 if ψ ˜ 1 ( t ) > 0 , U m a x if ψ ˜ 1 ( t ) < 0 ,
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 Section 3.2 and Section 3.3.

3.2. 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.
Using (16), (17), (19), and (20), we consider the following problem:
x ˙ ( t ) = A x ( t ) + B × max ( 0 , U m a x s i g n ( ψ 1 ( t ) ) ) , ψ ˙ ( t ) = A T ψ ( t ) , x ( 0 ) = ( 0 , 0 , 0 , 0 ) , x 1 ( t f ) = x e 1 , x 4 ( t f ) = x e 4 , ψ ( 0 ) is free , H ( t f , x ( t f ) , max ( 0 , U m a x s i g n ( ψ 1 ( t f ) ) ) , ψ ( t f ) ) = 0 .
Let z ( t ) = ( x ( t ) , ψ ( t ) ) . Then, we obtain the following two points’ boundary value problem:
z ˙ ( t ) = A * z ( t ) + B * , R ( z ( 0 ) , z ( t f ) ) = 0 ,
where A * M 8 × 8 ( R ) is the matrix given by
A * = A 0 4 × 4 0 4 × 4 A T ,
B * R 8 is the vector given by
B * = 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 if ψ 1 ( t ) > 0 , U m a x , 0 , 0 , 0 , 0 , 0 , 0 , 0 if ψ 1 ( t ) < 0 ,
and R ( z ( 0 ) , z ( t f ) ) is given by (2), (15), and (19). We consider the following Cauchy problem:
z ˙ ( t ) = A * z ( t ) + B * , z ( 0 ) = z 0 .
If we define the shooting function S : R 4 R 3 by
S ( z 0 ) = R ( t f , z ( t f , z 0 ) ) ,
where z ( t , z 0 ) is the solution of the Cauchy problem (25), then the two points’ boundary value problem (21) is equivalent to
S ( z 0 ) = 0 .
To solve (27), we use Newton’s method [23].

3.3. 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):
u ( t ) = U m a x , t [ 0 , t f ] .
Strategy 2 (zero switching times):
u ( t ) = 0 , t [ 0 , t f ] .
Strategy 3 (one switching time):
u ( t ) = U m a x if 0 t < t c , 0 if t c < t t f ,
where t c is a switching time.
Strategy 4 (one switching time):
u ( t ) = 0 if 0 t < t c , U m a x if t c < t t f .
Strategy 5 (two switching times):
u ( t ) = U m a x if 0 < t < t c 1 , 0 if t c 1 < t < t c 2 . U m a x if t c 2 < t t f ,
where t c 1 and t c 2 represent two switching times.
Strategy 6 (two switching times):
u ( t ) = 0 if 0 < t < t c 1 , U m a x if t c 1 < t < t c 2 . 0 if t c 2 < t t f .
Strategy 7 (three switching times):
u ( t ) = U m a x if 0 < t < t c 1 , 0 if t c 1 < t < t c 2 . U m a x if t c 2 < t t c 3 . 0 if t c 3 < t < t f ,
where t c 1 , t c 2 , and t c 3 represent three switching times.
Strategy 8 (three switching times):
u ( t ) = 0 if 0 < t < t c 1 , U m a x if t c 1 < t < t c 2 . 0 if t c 2 < t t c 3 . U m a x if t c 3 < t < t f .
Let x ( t ) be the trajectory associated with the control u ( t ) , given by the relation
x ( t ) = exp ( A t ) x ( 0 ) + 0 t exp ( A ( t s ) ) B u ( t ) d s ,
where exp ( A ) is the exponential matrix of A.
To calculate the switching times t c , t c 1 , t c 2 and the final time t f , we have to solve the following nonlinear equation:
x ˜ e F ( t f ) = ( x e 1 , x e 4 ) .
We also solve (37) using the Newton method [23].

4. 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.518 mg , 64.2371 mg , 813.008 mg , 3.4 mg ) , u e = 6.0907 mg / min .
Following the Schnider model, the matrix A of the dynamic system (11) is given by:
A = 0.9175 0.0683 0.0035 0 0.3020 0.0683 0 0 0.1960 0 0.0035 0 0.1068 0 0 0.4560 and B = 1 0 0 0 .
We are interested to solve the following minimum-time control problem:
min t f J = 0 t f d t , x ˙ ( t ) = A x ( t ) + B u ( t ) , x ( 0 ) = ( 0 , 0 , 0 , 0 ) , x e 1 ( t f ) = 14.518 mg , x e 4 ( t f ) = 3.4 mg , 0 u ( t ) 106.0907 , t [ 0 , t f ] , t f is free .

4.1. Numerical Resolution by the Shooting Method

Let z ( t ) = ( x ( t ) , ψ ( t ) ) . We consider the following Cauchy problem:
z ˙ ( t ) = A * z ( t ) + B * , z ( 0 ) = z 0 = ( 0 , 0 , 0 , 0 , ψ 01 , ψ 02 , ψ 03 , ψ 04 ) ,
where
A * = 10 4 9175 683 35 0 0 0 0 0 3020 683 0 0 0 0 0 0 196 0 35 0 0 0 0 0 1068 0 0 456 0 0 0 0 0 0 0 0 9175 3020 196 1068 0 0 0 0 683 683 0 0 0 0 0 0 35 0 35 0 0 0 0 0 0 0 0 456 ,
B * = max ( 0 , 106.0907 s i g n ( ψ 1 ( t ) ) ) 0 0 0 0 0 0 0 .
The shooting function S is given by
S ( z 0 ) = ( S 1 ( z 0 ) , S 2 ( z 0 ) , S 3 ( z 0 ) ) ,
where
S 1 ( z 0 ) = x e 1 ( t f ) 14.518 , S 2 ( z 0 ) = x e 4 ( t f ) 3.4 , S 3 ( z 0 ) = 1 + ψ T ( t f ) A x ( t f ) + B max ( 0 , 106.0907 s i n g ψ 1 ( t f ) ) .
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
t f = 1.8397 min ,
with
ψ T ( 0 ) = ( 0.0076 , 0.0031 , 0.0393 , 0.0374 ) .

4.2. 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 e F = ( 14.518 , 3.4 ) ,
that is,
u ( 0 ) = U m a x = 106.0907 mg / min .
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.0907 mg / min for all t [ 0 , t f ] . The trajectory x ( t ) , associated with this control u ( t ) , is given by the following relation:
x ( t ) = 0 t exp ( A ( t s ) ) B U m a x d s , t [ 0 , t f ] ,
where
exp ( A ( t s ) ) = V D ( t s ) V 1
with
V = 0 0.9085 0.0720 0.0058 0 0.3141 0.9377 0.0266 0 0.1898 0.3395 0.9996 1 0.1997 0.0187 0.0014
and
D ( τ ) = exp 0.4560 τ 0 0 0 0 exp 0.9419 τ 0 0 0 0 exp 0.0451 τ 0 0 0 0 exp 0.0024 τ .
System (37) takes the form
x 1 ( t f ) = 14.518 , x 4 ( t f ) = 3.4 ,
and has no solutions. Thus, Strategy 1 is not feasible.
Strategy 3: Let u ( t ) , t [ 0 , t f ] , be the control defined by
u ( t ) = 106.0907 mg / min if 0 t < t c , 0 if t c < t t f .
The trajectory x ( t ) associated with this control u ( t ) is given by
x ( t ) = 0 t exp ( A ( t s ) ) B U m a x d s if 0 t t c , exp ( A ( t t c ) ) x ( t c ) if t c < t t f ,
where
exp ( A ( t t c ) ) = V D ( t t c ) V 1 .
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
t c < t f .
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
t f = 1.8397 min ,
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
u ( t ) = 106.0907 mg / min if 0 t < t c 1 , 0 if t c 1 < t < t c 2 . 106.0907 mg / min if t c 2 < t t f ,
where t c 1 and t c 2 are the two switching times. The trajectory x ( t ) associated with control (59) is given by
x ( t ) = 0 t exp ( A ( t s ) ) B U m a x d s if 0 t t c 1 , exp ( A ( t t c 1 ) ) x ( t c 1 ) if t c 1 < t t c 2 , exp ( A ( t t c 2 ) ) x ( t c 2 ) + t c 2 t exp ( A ( t s ) ) B U m a x d s if t c 2 < t t f .
To compute the two switching times t c 1 and t c 2 and the final time t f , we have to solve the nonlinear system (52) with
0 t c 1 t c 2 t f .
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
u ( t ) = 106.0907 mg / min if 0 t < t c 1 , 0 if t c 1 < t < t c 2 . 106.0907 mg / min if t c 2 < t t c 3 , 0 mg / min if t c 3 < t t f ,
where t c 1 , t c 2 , and t c 3 are the three switching times. The trajectory x ( t ) associated with Control (62) is given by
x ( t ) = 0 t exp ( A ( t s ) ) B U m a x d s if 0 t t c 1 , exp ( A ( t t c 1 ) ) x ( t c 1 ) if t c 1 < t t c 2 , exp ( A ( t t c 2 ) ) x ( t c 2 ) + t c 2 t exp ( A ( t s ) ) B U m a x d s if t c 2 < t t c 3 , exp ( A ( t t c 3 ) ) x ( t c 3 ) if t c 3 < t t f .
To compute the three switching times t c 1 , t c 2 , and t c 3 and the final time t f , we have to solve the nonlinear system (52) with
0 t c 1 t c 2 t c 3 t f .
It turns out that System (52) subject to Condition (64) has no solution. Thus, Strategy 7 is also not feasible.
In Figure 2 and Figure 3, we present the solutions of the linear system of differential Equation (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 e F = ( 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.

5. 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 time-optimal 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.

Author Contributions

Conceptualization, M.A.Z., C.J.S., and D.F.M.T.; methodology, M.A.Z., C.J.S., and D.F.M.T.; software, M.A.Z.; validation, C.J.S. and D.F.M.T.; formal analysis, M.A.Z., C.J.S., and D.F.M.T.; investigation, M.A.Z., C.J.S., and D.F.M.T.; writing—original draft preparation, M.A.Z., C.J.S., and D.F.M.T.; writing—review and editing, M.A.Z., C.J.S. and D.F.M.T.; visualization, M.A.Z.; supervision, C.J.S. and D.F.M.T.; funding acquisition, M.A.Z., C.J.S., and D.F.M.T. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the Portuguese Foundation for Science and Technology (FCT—Fundação para a Ciência e a Tecnologia) through the R&D Unit CIDMA, Grant Numbers UIDB/04106/2020 and UIDP/04106/2020, and within the project “Mathematical Modelling of Multi-scale Control Systems: Applications to Human Diseases” (CoSysM3), Reference 2022.03091.PTDC, financially supported by national funds (OE) through FCT/MCTES.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

No new data were created or analyzed in this study. Data sharing is not applicable to this article. The numerical simulations of Section 4 were implemented in MATLAB R2022a. The computer code is available from the authors upon request.

Acknowledgments

The authors are grateful to four anonymous referees for their constructive remarks and questions that helped to improve the paper.

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

References

  1. Evers, A.S.; Maze, M.; Kharasch, E.D. Anesthetic Pharmacology: Basic Principles and Clinical Practice, 2nd ed.; Cambridge University Press: Cambridge, UK, 2011. [Google Scholar]
  2. Singh, N.; Sidawy, A.N.; Dezee, K.; Neville, R.F.; Weiswasser, J.; Arora, S.; Aidinian, G.; Abularrage, C.; Adams, E.; Khuri, S.; et al. The effects of the type of anesthesia on outcomes of lower extremity infrainguinal bypass. J. Vasc. Surg. 2006, 44, 964–970. [Google Scholar] [CrossRef]
  3. Merry, A.F.; Mitchell, S.J. Complications of anaesthesia. Anaesthesia 2018, 73, 7–11. [Google Scholar] [CrossRef]
  4. Beck, C.L. Modeling and control of pharmacodynamics. Eur. J. Control 2015, 24, 33–49. [Google Scholar] [CrossRef]
  5. Meibohm, B.; Derendorf, H. Basic concepts of pharmacokinetic/pharmacodynamic (PK/PD) modelling. Int. J. Clin. Pharmacol. Ther. 1997, 35, 401–413. [Google Scholar]
  6. Absalom, A.R.; Mani, V.; Smet, T.D.; Struys, M.M.R.F. Pharmacokinetic models for propofol—Defining and illuminating the devil in the detail. Br. J. Anaesth. 2009, 103, 26–37. [Google Scholar] [CrossRef] [PubMed]
  7. Enlund, M. TCI: Target Controlled Infusion, or Totally Confused Infusion? Call for an Optimised Population Based Pharmacokinetic Model for Propofol. Upsala J. Med. Sci. 2008, 113, 161–170. [Google Scholar] [CrossRef] [PubMed]
  8. Oshin, T.A. Exploratory mathematical frameworks and design of control systems for the automation of propofol anesthesia. Int. J. Dyn. Control 2022, 10, 1858–1875. [Google Scholar] [CrossRef]
  9. Wu, X.; Zhang, H.; Li, J. An analytical approach of one-compartmental pharmacokinetic models with sigmoidal Hill elimination. Bull. Math. Biol. 2022, 84, 117. [Google Scholar] [CrossRef]
  10. Nanditha, C.K.; Rajan, M.P. An adaptive pharmacokinetic optimal control approach in chemotherapy for heterogeneous tumor. J. Biol. Syst. 2022, 30, 529–551. [Google Scholar] [CrossRef]
  11. Wang, J.-J.; Dai, Z.; Zhang, W.; Shi, J.J. Operating room scheduling for non-operating room anesthesia with emergency uncertainty. Ann. Oper. Res. 2023, 321, 565–588. [Google Scholar] [CrossRef]
  12. Schnider, T.W.; Minto, C.F.; Gambus, P.L.; Andresen, C.; Goodale, D.B.; Shafer, S.L.; Youngs, E.J. The influence of method of administration and covariates on the pharmacokinetics of propofol in adult volunteers. Anesthesiology 1998, 88, 1170–1182. [Google Scholar] [CrossRef]
  13. Marsh, B.; White, M.; Morton, N.; Kenny, G.N. Pharmacokinetic model driven infusion of propofol in children. Br. J. Anaesth. 1991, 67, 41–48. [Google Scholar] [CrossRef] [PubMed]
  14. Zabi, S.; Queinnec, I.; Tarbouriech, S.; Garcia, G.; Mazerolles, M. New approach for the control of anesthesia based on dynamics decoupling. IFAC-PapersOnLine 2015, 48, 511–516. [Google Scholar] [CrossRef]
  15. Zabi, S.; Queinnec, I.; Garcia, G.; Mazerolles, M. Time-optimal control for the induction phase of anesthesia. IFAC-PapersOnLine 2017, 50, 12197–12202. [Google Scholar] [CrossRef]
  16. Ilyas, M.; Khan, A.; Khan, M.A.; Xie, W.; Riaz, R.A.; Khan, Y. Observer design estimating the propofol concentration in PKPD model with feedback control of anesthesia administration. Arch. Control Sci. 2022, 32, 85–103. [Google Scholar] [CrossRef]
  17. Bock, H.G.; Plitt, K.J. A Multiple shooting algorithm for direct solution of optimal control problems. IFAC Proc. Vol. 1984, 17, 1603–1608. [Google Scholar] [CrossRef]
  18. Pontryagin, L.S.; Boltyanskii, V.G.; Gamkrelidze, R.V.; Mishchenko, E.F. The Mathematical Theory of Optimal Processes; Trirogoff, K.N., Neustadt, L.W., Eds.; Interscience Publishers John Wiley and Sons, Inc.: New York, NY, USA, 1962. [Google Scholar]
  19. Bailey, J.M.; Haddad, W.M. Drug dosing control in clinical pharmacology. IEEE Control Syst. Mag. 2005, 25, 35–51. [Google Scholar] [CrossRef]
  20. James, W.P.T. Research on obesity: A report of the DHSS/MRC Group. Nutr. Bull. 1977, 4, 187–190. [Google Scholar] [CrossRef]
  21. Haddad, W.M.; Chellaboina, V.; Hui, Q. Nonnegative and Compartmental Dynamical Systems; Princeton Univ. Press: Princeton, NJ, USA, 2010. [Google Scholar]
  22. Zaitri, M.A.; Bibi, M.O.; Bentobache, M. A hybrid direction algorithm for solving optimal control problems. Cogent Math. Stat. 2019, 6, 1612614. [Google Scholar] [CrossRef]
  23. Deuflhard, P. Newton Methods for Nonlinear Problems; Springer Series in Computational Mathematics; Springer: Berlin, Germany, 2004; Volume 35. [Google Scholar]
  24. Elsisi, M.; Soliman, M. Optimal design of robust resilient automatic voltage regulators. ISA Trans. 2021, 108, 257–268. [Google Scholar] [CrossRef]
  25. Mohamed, M.A.E.; Mohamed, S.M.R.; Saied, E.M.M.; Elsisi, M.; Su, C.-L.; Hadi, H.A. Optimal energy management solutions using artificial intelligence techniques for photovoltaic empowered water desalination plants under cost function uncertainties. IEEE Access 2022, 10, 93646–93658. [Google Scholar] [CrossRef]
  26. Morales-Delgado, V.F.; Taneco-Hernández, M.A.; Vargas-De-León, C.; Gómez-Aguilar, J.F. Exact solutions to fractional pharmacokinetic models using multivariate Mittag-Leffler functions. Chaos Solitons Fractals 2023, 168, 113164. [Google Scholar] [CrossRef]
Figure 1. Schematic diagram of the PK/PD model with the effect site compartment of Bailey and Haddad [19].
Figure 1. Schematic diagram of the PK/PD model with the effect site compartment of Bailey and Haddad [19].
Axioms 12 00867 g001
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 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.
Axioms 12 00867 g002
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.
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.
Axioms 12 00867 g003
Figure 4. The optimal continuous infusion rate u ( t ) of the induction phase of anesthesia, as obtained by the shooting and analytical methods.
Figure 4. The optimal continuous infusion rate u ( t ) of the induction phase of anesthesia, as obtained by the shooting and analytical methods.
Axioms 12 00867 g004
Table 1. Parameter values for model (1) according to Schnider’s model [12].
Table 1. Parameter values for model (1) according to Schnider’s model [12].
ParameterEstimation
a 10 ( min 1 ) 0.443 + 0.0107 ( weight 77 ) 0.0159 ( LBM 59 ) + 0.0062 ( height 177 )
a 12 ( min 1 ) 0.302 0.0056 ( age 53 )
a 13 ( min 1 ) 0.196
a 21 ( min 1 ) 1.29 0.024 ( age 53 ) / 18.9 0.391 ( age 53 )
a 31 ( min 1 ) 0.0035
a e 0 ( min 1 ) 0.456
v 1 ( L ) 4.27
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Zaitri, M.A.; Silva, C.J.; Torres, D.F.M. An Analytic Method to Determine the Optimal Time for the Induction Phase of Anesthesia. Axioms 2023, 12, 867. https://doi.org/10.3390/axioms12090867

AMA Style

Zaitri MA, Silva CJ, Torres DFM. An Analytic Method to Determine the Optimal Time for the Induction Phase of Anesthesia. Axioms. 2023; 12(9):867. https://doi.org/10.3390/axioms12090867

Chicago/Turabian Style

Zaitri, Mohamed A., Cristiana J. Silva, and Delfim F. M. Torres. 2023. "An Analytic Method to Determine the Optimal Time for the Induction Phase of Anesthesia" Axioms 12, no. 9: 867. https://doi.org/10.3390/axioms12090867

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop