PRES22_0231.docx DOI: 10.3303/CET2294213 Paper Received: 16 April 2022; Revised: 09 May 2022; Accepted: 14 May 2022 Please cite this article as: Ai J., Dai J., Zhai C., Sun W., 2022, Pulse Drug Delivery Strategy Based on Pharmacokinetic-Pharmacodynamic (PK- PD) Model of Tolerance, Chemical Engineering Transactions, 94, 1279-1284 DOI:10.3303/CET2294213 CHEMICAL ENGINEERING TRANSACTIONS VOL. 94, 2022 A publication of The Italian Association of Chemical Engineering Online at www.cetjournal.it Guest Editors: Petar S. Varbanov, Yee Van Fan, Jiří J. Klemeš, Sandro Nižetić Copyright © 2022, AIDIC Servizi S.r.l. ISBN 978-88-95608-93-8; ISSN 2283-9216 Pulse Drug Delivery Strategy Based on Pharmacokinetic- Pharmacodynamic (PK-PD) Model of Tolerance Jiali Aia, Jindong Daia, Chi Zhaib,*, Wei Suna,* aBeijing University of Chemical Technology, Beijing, China bKunming University of Science and Technology, Kunming, China zhaichi@kmust.edu.cn and sunwei@ mail.buct.edu.cn In clinical practice, the curative effect of drugs will decrease due to the increase of drug resistance, especially in stable long-term drug delivery systems. Intermittent drug dosing is an effective method to reduce drug resistance and achieve a better effect comparing with consistent one. The relationship between drug concentration and effect can be described by pharmacokinetics-pharmacodynamics (PK-PD) model, and a more suitable drug delivery strategy can be obtained through analysis. In this paper, a periodic dosing strategy is proposed based on a simple PK-PD model of tolerance. Through the analysis of this model, a pulse drug delivery strategy is established, which is corresponding to the rhythm of human body and easy for drug design. By discretizing the differential equations into algebraic equations based on the collocation method, the drug delivery strategy is optimized. Then, the optimized drug delivery system is simulated by the fourth-order Runge-Kutta method. Based on this method, a reference for designing a drug delivery system can be provided to reduce unnecessary dosage and improve drug efficacy. 1. Introduction It is critical for clinic treatment to deliver drugs in a way that rationally matches the pharmacokinetic- pharmacodynamic (PK-PD) model. For many drugs, their effect depends not only on current concentration but also on exposure history. When the effect of the drug is weakened due to continuous exposure, patients are considered to have developed tolerance to the drug. Pochett et.al. (1988) proposed a PK-PD model that can be utilized to describe the development of tolerance. In this model, tolerance is attributed to the accumulation of drug metabolites as drug antagonists, which leads to the decrease of drug effects over time. On this basis, Varigonda et.al. (2004) proposed a dimensionless model in which all variables and drug inputs are non-negative. The (dimensionless) drug concentration and metabolite (antagonist) concentrations are denoted as 𝑐𝑐(0 ≤ c ≤ 10) and 𝑎𝑎(0 ≤ a ≤ 10). The pharmacokinetic model is shown in Eq(1). � 𝑑𝑑𝑑𝑑 𝑑𝑑𝑑𝑑 = 𝑢𝑢 − 𝑐𝑐 𝑑𝑑𝑑𝑑 𝑑𝑑𝑑𝑑 = 𝑘𝑘(𝑐𝑐 − 𝑎𝑎) (1) Where 𝑢𝑢(0 ≤ 𝑢𝑢 ≤ 10) represents the dosing rate, 𝑘𝑘 is a constant (usually 𝑘𝑘 = 0.1) that represents antagonist elimination rate. The drug effect is represented by 𝐸𝐸 and is given by Eq(2). 𝐸𝐸 = 𝑑𝑑 (1+𝑑𝑑)(1+𝑑𝑑 𝑑𝑑∗⁄ ) (2) There is an optimal stable solution which is given as 𝑢𝑢 = 𝑐𝑐 = 𝑎𝑎 = 𝑎𝑎∗ = 1. The optimal drug effect is given as 𝐸𝐸 = 0.25, as shown in Figure 1. However, according to pharmacodynamic, the drug's effective interval is 0.3 < 𝐸𝐸 < 0.6 in this model. Stable input cannot meet the therapeutic effect for a long time, which is unfavourable for some chronic diseases that require long-term medication. A classic example is nitroglycerin patch used to prevent angina pectoris. These patches should be worn 24 h/d and replaced with new patches daily. However, studies found that this protective effect diminished after about 12 h due to the development of drug resistance. The current dosing strategy 1279 recommends that the patient wear the patch for only 1/2 d and wait for resistance to recover (Parker et.al., 1995). This therapy increases exercise duration without significant evidence of nitrate tolerance or rebound phenomena for long time. Taking this case as an example, the simulation results of its drug concentration and drug effect are shown in the Figure 2, and it can provide effective drug effect within a certain period of time in each cycle. a b Figure 1: Optimal stable solution drug input strategy: (a) drug input rate and concentration changing, (b) drug effect. a b c Figure 2: Park input strategy: (a) drug input rate and concentration change, (b) drug effect, (c) concentration phase diagram. This is based on experimental observations, but if the PK-PD model is analysed, there are better drug delivery strategies. In order to find the optimal dosing trajectory, this problem is usually solved as an optimal period control problem with uncertain periods and waveforms. It is necessary to establish a performance index and then find out the maximum value of this index in this case. Varigonda et.al. (2004) defines such an index function in Eq(3), and the optimal object 𝐽𝐽 in Eq(4). They used Fourier basis functions to parameterize the input waveform. The optimal periodic dosing strategy was calculated using a flatness-based method and a shooting method. Compared with the current recommended periodic operation index about J = 0.3052, the optimized periodic operation index is about J = 0.3537, which improves the average therapeutic effect. Ghanaatpishe et.al. (2017) discussed the online periodic optimization method in the case of parameter uncertainty, and the optimized periodic operation index was also improved relative to the recommended periodic operation index. 𝐼𝐼 = (𝐸𝐸 𝐸𝐸1⁄ )𝛾𝛾 (1+(𝐸𝐸 𝐸𝐸1⁄ )𝛾𝛾)(1+(𝐸𝐸 𝐸𝐸2⁄ )2𝛾𝛾) (3) 𝐽𝐽 = 1 𝑇𝑇 ∫ 𝐼𝐼𝑇𝑇0 𝑑𝑑𝑑𝑑 (4) As an optimal control problem, these methods provide a better therapeutic effect than the periodic strategy observed in the experiment, but the rhythm of the human body is ignored in the established drug delivery system. The administration cycle that does not conform to the rhythm of the human body may lead to the change of pharmacodynamic relationship and the rejection of the patient. Moreover, this optimal dosing scheme obtained by the Fourier basis function has a large number of complex oscillations, which is difficult to achieve through drug design (Mclean and Zhan, 2021). This administration trajectory requires complex control equipment and cannot be achieved using currently commonly used sustained-release methods. Therefore, this problem is not an optimal period control problem with uncertain periods and waveform. In this paper, the problem of optimal periodic operation with piecewise function as the optimal input is presented. The paper is organized as follows: a pulse drug delivery system is introduced in the second section. The system 1280 is solved by collocation method and simulated by MATLAB in the third section. The fourth part is the conclusion of this paper. 2. Establishment of optimal pulse dosing strategy problems As a drug delivery system, its periodic changes should conform to the human rhythm. This means that when the system cycle is T = 24 h, it is most suitable. In order to facilitate drug design, the input drug concentration should be kept as stable as possible. Therefore, a drug delivery system with a 24 h period of pulsed drug trajectory will be introduced next. 2.1 Waveforms for periodic dosing For the system proposed by Eq(1) and Eq(2), steady dosing will lead to an increase in drug resistance and thus fail to achieve the minimum active drug effect. It is necessary to analyze the relation of 𝑐𝑐 and 𝑎𝑎 within the effective range of the drug, as shown in Figure 3. a b Figure 3. Relationship of 𝑐𝑐, 𝑎𝑎, 𝐸𝐸. a. Three-dimensional plot, b contour plot, where yellow represents the effective concentration range. It can be seen from Figure 3a that when antagonist concentration is high, the drug concentration has less influence on drug effect. This means that if the antagonist concentration is within a certain range, a particularly high drug concentration is required to produce an effective drug effect. The concentration of the antagonist should be controlled as low as possible from the viewpoint of reducing drug input. It can be seen from Figure 3b that when 𝑐𝑐 ∈ [3/7,10] and 𝑎𝑎 ∈ [0,67/33], the drug is in the effective range (yellow interval). Therefore, the concentration of the antagonist should not exceed the effective control range during the administration process Eq(5). 𝑐𝑐 0.6(1 + 𝑐𝑐) − 1 ≤ a ≤ 𝑐𝑐 0.3(1 + 𝑐𝑐) − 1 (5) According to Eq(1), 𝑎𝑎 tends to 𝑐𝑐 at a slower rate than 𝑐𝑐. In order to keep 𝑎𝑎 within the range, it is necessary to reduce 𝑐𝑐 to restore 𝑎𝑎 and make it at a smaller value, so that 𝐸𝐸 can reach the effective range within a period of time. In order to make 𝑎𝑎 decrease as fast as possible to enter the next period, 𝑐𝑐 should be decreased as soon as possible at this time. When 𝑢𝑢 = 0, 𝑐𝑐 decreases the fastest. Therefore, the waveform of the periodic dosing of the drug is as shown in Eq(6). 𝑢𝑢 = �𝑥𝑥, 0 < 𝑑𝑑 ≤ 𝜏𝜏 0, 24 ≥ 𝑑𝑑 > 𝜏𝜏 (6) This is a drug delivery system with a pulsatile dosing strategy. The system can be described as Eq(7) � 𝑑𝑑𝑐𝑐1 𝑑𝑑𝑑𝑑 = 𝑢𝑢 − 𝑐𝑐 𝑑𝑑𝑎𝑎1 𝑑𝑑𝑑𝑑 = 𝑘𝑘(𝑐𝑐 − 𝑎𝑎) ，0 < 𝑑𝑑 ≤ 𝜏𝜏；� 𝑑𝑑𝑐𝑐2 𝑑𝑑𝑑𝑑 = −𝑐𝑐 𝑑𝑑𝑎𝑎2 𝑑𝑑𝑑𝑑 = 𝑘𝑘(𝑐𝑐 − 𝑎𝑎) ，𝜏𝜏 < 𝑑𝑑 ≤ 24 (7) 2.2 Establishment of periodic system of drug delivery The system needs to meet the effective requirements, that means 𝐸𝐸𝑚𝑚𝑑𝑑𝑚𝑚 > 0.3. According to Eq(6), 𝐸𝐸(𝑑𝑑) = 0.3 should only have two solutions in a period. 𝛼𝛼 and 𝛽𝛽 should satisfy Eq(9). The certification process is in the appendix. E(α) = E(𝛽𝛽) = 0.3，∆𝑑𝑑 = 𝛽𝛽 − α (8) α < β and β = 𝜏𝜏 (9) Although 𝑢𝑢 is a piecewise function, 𝑎𝑎 and 𝑐𝑐 are continuous functions, there should be Eq(10) 1281 𝑐𝑐1(𝜏𝜏) = 𝑐𝑐2(𝜏𝜏); 𝑎𝑎1(𝜏𝜏) = 𝑎𝑎2(𝜏𝜏) (10) According to Eq(9), at time 𝜏𝜏, 𝑎𝑎 and 𝑐𝑐 should meet Eq(11) 𝑎𝑎(𝜏𝜏) = 𝑐𝑐(𝜏𝜏) 0.3�1 + 𝑐𝑐(𝜏𝜏)� − 1 (11) Although the performance index Eq (3) proposed by Varigonda et.al. (2004) provides a method to solve the optimal cycle control under the condition of cycle uncertainty, and has achieved appropriate results, it is not the ideal optimal parameter in practice. ∆𝑑𝑑 is chosen as the optimal parameter, the periodic optimization problem can be described as Eq(12). 𝑚𝑚𝑎𝑎𝑥𝑥 𝐽𝐽 =∆𝑑𝑑 𝑠𝑠. 𝑑𝑑. 𝑑𝑑𝑐𝑐1 𝑑𝑑𝑑𝑑 = 𝑢𝑢 − 𝑐𝑐; 𝑑𝑑𝑎𝑎1 𝑑𝑑𝑑𝑑 = 𝑘𝑘(𝑐𝑐 − 𝑎𝑎); 𝑐𝑐(0) = 𝑐𝑐(𝑇𝑇);𝑎𝑎(0) = 𝑎𝑎(𝑇𝑇) (12) 3. Numerical result for periodic dosing strategy The optimal periodic control problem needs to be determined by the π-test. However, the π-test considers the case of a continuous function, and this part needs to find the best value of the pulse period problem. Therefore, it can not directly use the optimal periodic control method to solve this problem, and the periodic dosing strategy will be discussed below through the collocation method. 3.1 Theory of calculating periodic solutions by collocation method It is assumed that the general autonomous system has periodic solutions with period 𝑇𝑇, and its mathematical problem is shown in Eq(13). 𝑑𝑑𝑥𝑥 𝑑𝑑𝑑𝑑 = 𝑓𝑓(𝑥𝑥) 𝑥𝑥|𝑑𝑑=0 = 𝑥𝑥|𝑑𝑑=𝑇𝑇 (13) where 𝑥𝑥 ∈ 𝑅𝑅𝑛𝑛 is an N-dimensional variable matrix, 𝑓𝑓(𝑥𝑥): 𝑅𝑅𝑛𝑛 → 𝑅𝑅𝑛𝑛 is a nonlinear mapping. It can become a boundary value problem (BVP) as Eq(14). � 𝑑𝑑𝑥𝑥 𝑑𝑑𝜏𝜏 = 𝑇𝑇𝑓𝑓(𝑥𝑥) 𝑥𝑥|𝜏𝜏=0 = 𝑥𝑥|𝜏𝜏=1, 𝜏𝜏 = 𝑑𝑑/𝑇𝑇 ∈ [0, 1] (14) Due to the piecewise input function, the model cannot obtain a numerically unique solution. Shooting method and Finite difference method are the common methods for calculating BVP problem. But considering the periodic solution of initial value sensitivity, this paper discusses the application of Collocation method. The idea is to approximate the ode in Eq(13) at a specific collocation point, thus discretizing the original problem. For the BVP model of a given autonomous system, the optimal approximation method is used to obtain the optimal numerical approximation under the given calculation amount. Figure 4. Location intervals and collocation points As shown in Figure 4, assuming at a ≤ t ≤ b, there are n collocation points, in which the 𝑖𝑖 discrete interval [𝑑𝑑𝑖𝑖 , 𝑑𝑑𝑖𝑖 + 1 ] with 𝑘𝑘 − 1 order polynomial 𝑥𝑥𝜋𝜋 interpolation conditions (Eq(15)). 𝑥𝑥𝜋𝜋(𝑑𝑑𝑖𝑖) = 𝑥𝑥𝑖𝑖 , 𝑥𝑥𝜋𝜋′(𝑑𝑑𝑖𝑖) = 𝑓𝑓�𝑥𝑥𝑖𝑖,𝑗𝑗 ,𝑢𝑢�, 1 ≤ 𝑖𝑖 ≤ 𝑚𝑚, 1 ≤ 𝑗𝑗 ≤ 𝑛𝑛 (15) 𝑥𝑥𝜋𝜋(𝑑𝑑) approximate solution can be through the initial value for 𝑥𝑥𝜋𝜋(𝑑𝑑𝑖𝑖) numerical description (Eq(16)). 𝑥𝑥𝜋𝜋(𝑑𝑑) = 𝑥𝑥𝜋𝜋(𝑑𝑑𝑖𝑖) + ∫ ∑ 𝑓𝑓�𝑥𝑥𝑖𝑖,𝑗𝑗 ,𝑢𝑢𝑖𝑖,𝑗𝑗�𝑛𝑛 𝑗𝑗=1 𝐿𝐿𝑗𝑗(𝜉𝜉)𝑑𝑑 𝑑𝑑𝑖𝑖 𝑑𝑑𝜉𝜉 (16) where 𝑥𝑥𝑖𝑖,𝑗𝑗 = 𝑥𝑥𝜋𝜋�𝑑𝑑𝑖𝑖，𝑗𝑗�, 𝐿𝐿𝑗𝑗(𝜉𝜉) = ∏ 𝜉𝜉−𝜉𝜉𝑘𝑘 𝜉𝜉𝑗𝑗−𝜉𝜉𝑘𝑘 𝑛𝑛 𝑘𝑘=1,𝑘𝑘≠𝑗𝑗 , 𝜉𝜉 = 𝑑𝑑−𝑑𝑑𝑖𝑖 ℎ𝑖𝑖 , 𝜉𝜉 ∈ [0,1], 𝜉𝜉𝑗𝑗 = 𝑑𝑑𝑖𝑖,𝑗𝑗−𝑑𝑑𝑖𝑖 ℎ𝑖𝑖 , ℎ𝑖𝑖 = 𝑑𝑑𝑖𝑖+1 − 𝑑𝑑𝑖𝑖 (17) The Gaussian integral is taken as the integral of [0, 1], and the integral point is taken as the collocation point. Two sets of discrete equations (Eq(18)) are obtained. 𝑥𝑥𝜋𝜋�𝑑𝑑𝑖𝑖,𝑗𝑗� = 𝑥𝑥𝑖𝑖 + ℎ𝑖𝑖 ∑ 𝛼𝛼𝑗𝑗𝑘𝑘𝑓𝑓(𝑥𝑥𝑖𝑖𝑘𝑘)𝑛𝑛 𝑘𝑘=1 (18) 𝑥𝑥𝜋𝜋(𝑑𝑑𝑖𝑖+1) = 𝑥𝑥𝑖𝑖 + ℎ𝑖𝑖 ∑ 𝛽𝛽𝑗𝑗𝑘𝑘𝑓𝑓(𝑥𝑥𝑖𝑖𝑘𝑘)𝑛𝑛 𝑘𝑘=1 where 1282 𝛼𝛼𝑗𝑗𝑘𝑘 = ∫ 𝐿𝐿𝑗𝑗(𝜉𝜉)𝑑𝑑𝜉𝜉𝜉𝜉𝑘𝑘 0 ,𝛽𝛽𝑗𝑗𝑘𝑘 = ∫ 𝐿𝐿𝑗𝑗(𝜉𝜉)𝑑𝑑𝜉𝜉1 0 (19) By choosing the orthogonal configuration, Eq(18) can have the same matrix form as Runge-Kutta algorithm. The 𝛼𝛼𝑗𝑗𝑘𝑘 and 𝛽𝛽𝑗𝑗𝑘𝑘 in Eq(21) and the collocation point 𝜉𝜉 can be calculated by using Eq(19). To compute periodic solutions for a wide range of parameter values, the largest possible step size is employed. The simulation method is fourth-order Runge-Kutta method. It is given by Eq(20). 𝑥𝑥𝑛𝑛+1 = 𝑥𝑥𝑛𝑛 + ℎ 6 (𝑘𝑘1 + 2𝑘𝑘2 + 2𝑘𝑘3 + 𝑘𝑘4) (20) where 𝑘𝑘1 = 𝑓𝑓(𝑑𝑑𝑛𝑛, 𝑥𝑥𝑛𝑛), 𝑘𝑘2 = 𝑓𝑓 �𝑑𝑑𝑛𝑛 + ℎ 2 , 𝑥𝑥𝑛𝑛 + ℎ 2 𝑘𝑘1� , 𝑘𝑘3 = 𝑓𝑓 �𝑑𝑑𝑛𝑛 + ℎ 2 , 𝑥𝑥𝑛𝑛 + ℎ 2 𝑘𝑘2� , 𝑘𝑘4 = 𝑓𝑓(𝑑𝑑𝑛𝑛 + ℎ, 𝑥𝑥𝑛𝑛 + ℎ𝑘𝑘3) (21) 3.2 Numerical simulation results of impulsive periodic dosing strategy According to the collocation method, the model can be obtained to obtain a pulsatile periodic drug delivery strategy - Eq(22). 𝑢𝑢 = �1.2, 0 < 𝑑𝑑 ≤ 10 0, 24 ≥ 𝑑𝑑 > 10 (22) The numerical simulation results are shown in Figure 5. The 𝑚𝑚𝑎𝑎𝑥𝑥 𝐽𝐽 =∆𝑑𝑑 = 9.43. a b c d Figure 5: Impulsive periodic dosing strateg: (a) drug infusion rate (b) 𝑐𝑐 and 𝑎𝑎 concentration changge, (c) drug effect, (d) concentration phase diagram. As can be seen in Figure 5, the effective concentration of the drug and the concentration of the antagonist form a limit cycle. This shows that the scheme has formed a stable cycle. Meanwhile, in each cycle, there is a 9.43 h onset time. The results of steady-state input are not discussed here, as it is evident that there is no long-term stable periodic drug effect. Compared with the experimental experience method mentioned in the introduction (Figure 2), the 𝑚𝑚𝑎𝑎𝑥𝑥 𝐽𝐽 =∆𝑑𝑑 = 8.25. After analysis and optimization, the effective action time of the drug in each cycle (one day) was prolonged by 1.18 h. Compared with the optimized results of other literature, the average curative effect of each cycle is relatively low. However, from the effective time of each cycle, the effective time of other literature is only 6 - 8 h per cycle, which is lower than the optimization result. The reason for this result is that other literature shorten the time of each cycle, so as to improve the average curative effect of the cycle. Because when the drug input rate is greater than 1, the drug effect will reach its maximum value within 6 h. Continuous infusion will only lead to a rapid increase in antagonist concentration, thus reducing the curative effect. Their period is usually 15 - 18 h. This kind of processing is suitable for controlling the results. But for patients, the rapid cycles turnover means that more operations are needed. For systems that require long-term medication, patients may have a rebellious attitude and refuse to treat them. In addition, no matter what the waveform of the drug input is, there is no possibility of long-term stable and effective drug effect. This means that other drugs need to make up for the lack of efficacy during the period when the drug is not working. The segmented input provided herein is more suitable for co-administration. 1283 4. Conclusions In this paper, a pulsatile drug delivery system is proposed based on a PK-PD model of tolerance. This is not a cycle solution with the best average efficacy, but a dosing regimen that conforms to the rhythm of the body. At the same time, this pulsatile dosing regimen is also easy to achieve through the design of sustained drug delivery. The specific optimal input result obtained by the collocation method is shown in Eq(23). Stable periodic results are obtained based on the model. The results of this study provide an easy-to-realize solution to solve the decrease of clinic effectiveness due to drug-resistance, therefore keeping relatively low dosage and control potential side effect to patient. Less drug demand will keep the capacity of pharmaceutical production and related emissions. It is also a contribution to sustainable development of our society. Appendix During the time when 𝑢𝑢 is steadily inputting 𝜏𝜏 with quantity 𝑥𝑥, 𝐸𝐸 = 𝑒𝑒𝑘𝑘𝑑𝑑(−1 + 𝑘𝑘)((−1 + 𝑒𝑒𝑑𝑑)𝑥𝑥 + 𝑐𝑐(0)) (−𝑥𝑥 + 𝑒𝑒𝑑𝑑(1 + 𝑥𝑥) + 𝑐𝑐(0))(𝑥𝑥 + 𝑒𝑒𝑘𝑘𝑑𝑑(−1 + 𝑘𝑘)(1 + 𝑥𝑥) − 𝑒𝑒(−1+𝑘𝑘)𝑑𝑑𝑘𝑘(𝑥𝑥 − 𝑐𝑐(0)) − 𝑘𝑘𝑐𝑐(0)) (23) 𝑑𝑑𝐸𝐸 𝑑𝑑𝑑𝑑 = 𝑒𝑒𝑘𝑘𝑑𝑑(−1 + 𝑘𝑘) 𝐵𝐵(𝑑𝑑) − 𝐶𝐶(𝑑𝑑) 𝐴𝐴2 (24) where 𝐴𝐴 = (−𝑥𝑥 + 𝑒𝑒𝑑𝑑(1 + 𝑥𝑥) + 𝑐𝑐(0))2(𝑥𝑥 + 𝑒𝑒𝑘𝑘𝑑𝑑(−1 + 𝑘𝑘)(1 + 𝑥𝑥) − 𝑒𝑒(−1+𝑘𝑘)𝑑𝑑𝑘𝑘(𝑥𝑥 − 𝑐𝑐(0)) − 𝑘𝑘𝑐𝑐(0))2 𝐵𝐵(𝑑𝑑) = 2𝑒𝑒𝑘𝑘𝑑𝑑𝑘𝑘𝑥𝑥(𝑥𝑥 − 𝑐𝑐(0))2 + 𝑘𝑘(𝑥𝑥 − 𝑐𝑐(0))2(𝑥𝑥 − 𝑘𝑘𝑐𝑐(0)) + 𝑒𝑒2𝑑𝑑𝑘𝑘𝑥𝑥(1 + 𝑥𝑥)(𝑥𝑥 − 𝑘𝑘𝑐𝑐(0)) + 𝑒𝑒𝑑𝑑(1 − 𝑘𝑘 − 2𝑘𝑘𝑥𝑥)(𝑥𝑥 − 𝑐𝑐(0))(𝑥𝑥 − 𝑘𝑘𝑐𝑐(0)) 𝐶𝐶(𝑑𝑑) = 𝑒𝑒𝑑𝑑+𝑘𝑘𝑑𝑑(1 + 𝑘𝑘(−1 + 𝑥𝑥))(1 + 𝑥𝑥)(𝑥𝑥 − 𝑐𝑐(0)) + 𝑒𝑒(−1+𝑘𝑘)𝑑𝑑𝑘𝑘(𝑥𝑥 − 𝑐𝑐(0))3 (25) It can be seen that when 𝐵𝐵(𝑑𝑑) > 𝐶𝐶(𝑑𝑑), E is a monotonically increasing function. When 𝐴𝐴 ≠ 0 and 𝐵𝐵(𝑑𝑑) = 𝐶𝐶(𝑑𝑑) , 𝐸𝐸 takes the most extreme value. The problem becomes a problem analyzing 𝐵𝐵(𝑑𝑑) and 𝐶𝐶(𝑑𝑑). As shown in Figure 6, when 𝑥𝑥 ∈ [1,2], 𝑐𝑐(0) ∈ [0,1], 𝐵𝐵(𝑑𝑑) and 𝐶𝐶(𝑑𝑑) is a monotonically increasing function and E has one and only one extreme point. a b c Figure 6: Schematic diagram of the change of 𝐵𝐵(𝑑𝑑) and 𝐶𝐶(𝑑𝑑) with time. u=1, 𝑐𝑐(0) = 𝑎𝑎(0) = 0 To maximize the optimization objective ∆𝑑𝑑, 𝑚𝑚𝑎𝑎𝑥𝑥 𝐸𝐸 must be greater than 0.3. This means that E(𝛽𝛽) = 0.3 is in the monotonically decreasing part. The continued input of 𝑢𝑢 will decrease drug effect due to the continuous increase of 𝑎𝑎. It is better to stop the input and let 𝑎𝑎 to recover. Acknowledgments The National Natural Science Foundation of China (Grant No. 21878012). References Ghanaatpishe M., Kehs M.A., Fathy H.K., 2017, Online Shaping of a Drug's Periodic Administration Trajectory for Efficacy Maximization, IEEE Transactions on Control Systems Technology, 25(6),1989-1999. Mclean K., Zhan W., 2021, Mathematical modelling of nanoparticle-mediated topical drug delivery to skin tissue, International Journal of Pharmaceutics, 611, 121322 Parker J.O., Amies M.H.,Hawkinson R.W., Heilman J.M., Hougham A.J., Vollmer M.C., Wilson R.R., 1995, Intermittent transdermal nitroglycerin therapy in angina pectoris, clinically effective without tolerance or rebound, Circulation, 91(5), 1368-1374 Porchet H.C. , Benowitz N.L., Sheiner L.B., 1988, Pharmacodynamic model of tolerance: Application to nicotine, Journal of Pharmacology and Experimental Therapeutics, 244(1), 231-236. Varigonda S. , Georgiou T.T., Siegel R.A., Daoutidis P., 2008, Optimal periodic control of a drug delivery system, Computers & Chemical Engineering, 32(10), 2256-2262. 1284 PRES22_0491.pdf Pulse Drug Delivery Strategy Based on Pharmacokinetic-Pharmacodynamic (PK-PD) Model of Tolerance