Modeling, Parameter Identification, and Adaptive Control of

a major contributor to the number one (acute myocardial infarc- tion), number ... timizing the administration of the anticoagulant drug heparin. The a...
1 downloads 0 Views 889KB Size
24 Modeling, Parameter Identification, and Adaptive Control of Anticoagulant Drug Therapy

Downloaded by EAST CAROLINA UNIV on January 2, 2018 | http://pubs.acs.org Publication Date: May 30, 1980 | doi: 10.1021/bk-1980-0124.ch024

THOMAS J. McAVOY Department of Chemical Engineering, Goessman Laboratory, University of Massachusetts, Amherst, MA 01003

The scope of the medical problems which require the use of anticoagulant drugs cannot be overstated. Arterial thrombosis is a major contributor to the number one (acute myocardial infarction), number three (stroke) and number four (renal) causes of deathinthe United States (1). Venous thromboembolismisthe most common non-surgical cause of deathinpatients hospitalized for major orthopedic procedures, the most frequent non-obstetrical cause of postpartum death, and a major cause of deathinpatients with chronic cardiac and pulmonary disease (1). Venous thrombosis is estimated to lead to the hospitalization of approximately 300,000 patients annuallyinthe United States, of which more than 50,000 die (1). The average hospital stay for the 250,000 survivors is a week to 10 days and the cost of this hospitalization has been estimated as $750 million (1). If a one day reduction in hospital stay could be achieved through optimumdosingof anticoagulants, then the potential savings per year can be estimated as $50 million. It has further been estimated (2) that through improved administration of anticoagulants 8,000 to 12,000 lives could be saved each year in the United States alone. The purpose of this paperisto illustrate how mathematical modeling and control theory can be applied to the problem of optimizing the administration of the anticoagulant drug heparin. The approach presented has the potential of solving some of the economic and life saving problems discussed above. Mechanism of Clot Formation and Drug Action. Blood clotting is an extremely complicated process which involves an enzyme cascade. To date 13 enzymes, called clotting factors, have been identified and these are normally presentinbloodintheir inactive state. Through an injury of some type one of the early enzymesisactivated. This enzymeinturn catalyses the activation of the next enzymeinthe cascade. Throughout the process a large amplification occurs. One molecule of a clotting factoriscapable of activating many molecules of the next enzymeinthe sequence and so forth. The final step in blood clotting involves 0-8412-0549-3/80/47-124-425$05.00/0 © 1980 American Chemical Society Squires and Reklaitis; Computer Applications to Chemical Engineering ACS Symposium Series; American Chemical Society: Washington, DC, 1980.

Downloaded by EAST CAROLINA UNIV on January 2, 2018 | http://pubs.acs.org Publication Date: May 30, 1980 | doi: 10.1021/bk-1980-0124.ch024

426

C O M P U T E R APPLICATIONS T O C H E M I C A L ENGINEERING

the polymerization of fibrinogen to fibrin. This reaction is catalysed by the clotting factor thrombin. Obviously the ability of blood to clot is extremely important when one has a wound or cut. However when clots form inside a vein or an artery a serious problem exists. Anticoagulant drugs are given to prevent this latter situation. Recent experimental evidence (3^4.) has demonstrated that heparin's anticoagulant effect results from its ability to catalyse the deactivation of several of the clotting factors. Of particular importance and the primary reason for heparin's anticoagulant effect is the drug's ability to deactivate thrombin by catalysing the reaction between thrombin and its inhibitor antithrombin III (AT-III). Because it acts directly to deactivate clotting factors, heparin has a very fast anticoagulant effect. Laboratory Measurement and Goals of Anticoagulant Therapy. For in-hospital patients receiving heparin, a blood sample is drawn, at least daily. The red cells are centrifuged from the sample leaving only plasma. An agent that activates the clotting cascade is added to a fixed amount of plasma. The time required for a clot to form, called the clotting time, is then measured. In addition to the patient's plasma, normal plasma is also tested. The goal of anticoagulant therapy is to have the patient's clotting time be 2 to 2-1/2 times that for the normal plasma. For heparin the most common clotting test is called the partial thromboplastin time, PTT, Modeling and Optimal Control of Heparin Modeling Heparin Itself. For heparin it is reasonable to use the CSTR model shown in Figure 1, Drug can be inputted to the reactor either by a bolus injection (impulse forcing) or by continuous intravenous infusion, Within the reactor drug is continuously removed from the blood, However, the mechanism of this removal is at present unknown. As discussed by McAvoy (5) there are some interesting dynamic properties associated with "the response of hep^ arin itself. Figure 2 shows a plot of experimental data on heparin together with the response of one model that has been proposed for the drug (j>). To generate the experimental data, three different bolus injections of heparin were given to healthy volunteers. Blood samples were then drawn periodically and heparin concentrations measured. When the logarithm of heparin concentration is plotted versus time, an apparently linear response is obtained. However, the slope of the response changes with each bolus. This change in slope indicates that an unusual nonlinear phenomena takes place when humans eliminate heparin. Stated mathematically, the initial condition inside the CSTR somehow is stored and it affects the entire transient response of the drug. The model response plotted in Figure 2 is based on the following metabolite-inhibition mechanism {§)

Squires and Reklaitis; Computer Applications to Chemical Engineering ACS Symposium Series; American Chemical Society: Washington, DC, 1980.

Downloaded by EAST CAROLINA UNIV on January 2, 2018 | http://pubs.acs.org Publication Date: May 30, 1980 | doi: 10.1021/bk-1980-0124.ch024

Figure 1.

CSTR model for heparin

Squires and Reklaitis; Computer Applications to Chemical Engineering ACS Symposium Series; American Chemical Society: Washington, DC, 1980.

COMPUTER

Downloaded by EAST CAROLINA UNIV on January 2, 2018 | http://pubs.acs.org Publication Date: May 30, 1980 | doi: 10.1021/bk-1980-0124.ch024

428

APPLICATIONS

T O CHEMICAL

ENGINEERING

.1.08-

.05·

.02-

.01 ~I20

îeÔ

240

300

360

420

TIME-MINUTES

Figure 2. Heparin data and metabolite-inhibition model for three separate injections of drug. ( ) heparin concentration; (O) experimental data; ( ) metabolite concentration, (a) 400 μ/kg bolus; (b) 200 μ/kg bolus; (c) 100 μ,/kg bolus (5).

Squires and Reklaitis; Computer Applications to Chemical Engineering ACS Symposium Series; American Chemical Society: Washington, DC, 1980.

24.

Anticoagulant

MCAVOY

Drug

Therapy

429

H + E £ HE + M-j+E M

1

+ Ε t M

Downloaded by EAST CAROLINA UNIV on January 2, 2018 | http://pubs.acs.org Publication Date: May 30, 1980 | doi: 10.1021/bk-1980-0124.ch024

1

(1 )

M-|E

(2)

- M

(3)

2

Heparin, H, is assumed to be metabolized by enzyme Ε to product M-. which competes with the drug for binding sites on E, In addition, M-. is biotransformed to which is non-competitive. If the mech­ anism given by equations f-3 is assumed and balances on Η and M, are made for the reactor shown in Figure 1, the following model for heparin results (6j 3Έ

I-V.H/iH + K j l + i y K p ) )

=

(4)

dM, +

V

1+

K

dTV P» - ep 1 where I is the rate of drug infusion divided by the reactor volume. Equations 4 and 5 contain four constants, V , Κ , K and k , These constants were fitted to the experimental data shown τη Figure 2 by numerically integrating the differential equations and then using a least square simplex approach on the results (7_). The values obtained are: V = .938, Κ = 15.7, K = ,385, and k = .00187. The responses predicted by equations 4 and 5 are alio shown in Figure 2. Note that 1=0 for these responses. As can be seen, equations 4 and 5 do an excellent job of fitting the ex­ perimental data. The M, response is given as the dashed line in Figure 2. As can be seen, M, builds up very quickly and achieves an almost constant value for most of the Η transient. In fact it is M, which "remembers" the initial condition in the reactor. In effedt M-| blocks an approximately constant number of binding sites on Ε for each response and this produces the apparent first order response. Since the number of blocked sites increases with the initial condition for H, the slope of each of the Η responses is different. In addition to the metabolite-inhibition model, McAvoy also presents an alternative model for heparin involving two non­ linear differential equations and based on a phagocytosis elimina­ tion mechanism (5J. =

k

M

5

p

p

Modeling the Clotting Time. Since the main objective of hep­ arin therapy is to control a patient's clotting time (normally the PTT), it is necessary to have a model for the clotting time. Both the nature of the clotting cascade and published experimental re­ sults (8) indicate that clotting time is not a simple function of heparin concentration alone. Other variables which affect the clotting time include the concentration and degree of activation of the various clotting factors and the concentration of other pro­ teins in the blood, such as antithrombin III. There is a degree

Squires and Reklaitis; Computer Applications to Chemical Engineering ACS Symposium Series; American Chemical Society: Washington, DC, 1980.

Downloaded by EAST CAROLINA UNIV on January 2, 2018 | http://pubs.acs.org Publication Date: May 30, 1980 | doi: 10.1021/bk-1980-0124.ch024

430

COMPUTER

APPLICATIONS

TO

CHEMICAL

ENGINEERING

of uncertainty at present as to what models will ultimately prove useful in describing heparin's effect on the clotting time because the necessary experimental data have not yet been gathered. The author is involved in a research program to generate these data. To illustrate the type of modeling approach that will be studied as well as the control problem which results, one specific example involving antithrombin III will be discussed in detail. This ex­ ample is typical of the modeling which will be carried out. Rosenberg's data (4) indicate that heparin binds to AT-III, while those of Marciniak (9) show that such binding does not take place. Both authors, however, agree that heparin's main route of action is via the enhancement of the inhibitory reaction which takes place between AT-III and thrombin. Marciniak and Gockerman (10) have also recently published data on AT-III for 24 patients receiving heparin infusions. These authors found that heparin therapy produced an exponential like transient decay in AT-III levels. After 3 days of heparin therapy AT-III levels approached a constant value which was 30% lower than normal. After cessation of heparin, AT-III levels returned to normal in about 3 days. Since AT-III is so important to heparin's action, it is pro­ posed to model its dynamic behavior, and to attempt to incorporate its concentration into an equation for clotting time. Thus, one equation which is proposed for describing clotting time is Ε = g (H, AT-III)

(6)

where Ε, the relative extension of the clotting time, is defined as Ε = (PTT/PTT - 1)

(7)

n

and PTT is the clotting time and PTT is the normal clotting time. The reason for considering an algebraic relationship between E, AT-III,and.Η is that the results of Rosenberg (4) indicate that heparin almost instantaneously catalyses the deactivation of throm­ bin by AT-III. Free thrombin levels in turn affect E. Thus, changes in both heparin level and AT-III level should have a very rapid effect on E. In order to make equation 6 more specific it will be assumed that heparin binds to antithrombin III n

Η + AT-III t

(H-AT-III)

(8)

Further, it will be assumed that Ε is proportional to the concen­ tration of the heparin antithrombin III cofactor Ε = m (H-AT-III)

(9)

If thermodynamic equilibrium is assumed for equation 8, equation 9 can be written as

Squires and Reklaitis; Computer Applications to Chemical Engineering ACS Symposium Series; American Chemical Society: Washington, DC, 1980.

24.

MCAVOY

Anticoagulant

Ε = m K

Drug

Therapy

431

(H)(AT-III)

eq

(10)

where Κ is the equilibrium constant for equation 8, To complete the modër it is necessary to have an equation for antithrombin III. One model which Marciniak and Gockerman's data (10) suggest will be appropriate is (AT-III) = - k (AT-III - AT-III ) - k H

Downloaded by EAST CAROLINA UNIV on January 2, 2018 | http://pubs.acs.org Publication Date: May 30, 1980 | doi: 10.1021/bk-1980-0124.ch024

ft

0

g

(11)

where k. and kg are constants and AT-III is the normal blood concentration of antithrombin III. At steaay state with no heparin present equation 11 forces AT-III to return to its steady state value. From Marciniak and Gockerman's data a value of AT-III can be gotten as 36 mg/d£. To achieve a 3 day transient response for AT-III a reasonable value of k. to use is .001 min. Lastly, if a steady state heparin concentration of .5 units/ml plasma produces the 30% reduction in AT-III levels then k can be calculated as .0216 ((mg/d*)/(units/ml min)) by setting d(AT-III)/dt = 0. R D

Optimal Control of Heparin Therapy. One complete model for a patient's response to heparin is given by equations 4, 5, 10 and 11. All of the parameters in the model have been specified except m Κ in equation 10. A reasonable value for m Κ can be gotten byassuming that Ε is 1,0 when Η is .5 and AT-IIT is 70% of AT-III . This assumption gives m K as .0794. Using this model and modern control theory it is possible to formulate an optimum dosing policy. Prior to heparin being administered the initial conditions in the CSTR are that AT-III equals AT-III and both Η and M, are 0. One clinically reasonable performance index is to choose the heparin infusion, I, to minimize the time that it takes to get from this initial state to an anticoagulated state where the clotting time is twice normal (E = 1). Mathematically one should minimize the performance index q

P.I. = J

e f Q

dt

(12)

subject to E(e ) = 1 = .0794 (H)(AT-III) f

(13)

In addition, I would be constrained to be less than some maximum clinically acceptable value I* I < I*

(14)

A reasonable value for I* is .968 units/hr ml plasma. This value results from a maximum infusion rate of 3000 units/hr into a typ­ ical plasma volume of 3100 ml. The problem posed involves time optimal control to a region. The approach which can be used to

Squires and Reklaitis; Computer Applications to Chemical Engineering ACS Symposium Series; American Chemical Society: Washington, DC, 1980.

432

COMPUTER

TO

CHEMICAL

ENGINEERING

solve this problem has been discussed by Denn (11). By appending 3 adjoint equations to the state equations, one can start at the final time, θ^, and numerically integrate backwards in time to solve for the optimal policy. To carry out the integration it is necessary to estimate H(e ) and M-.(e) and iterate on these values until the initial conditions on M; and AT-III are satisfied. The time optimal policy calls for an initial injection of 1420 units of heparin with no infusion (I = 0) for 5.25 minutes. Once a pa­ tient is brought to an anticoagulated state it is necessary to keep him there and still not violate equation 14, McAvoy {5) has presented an infusion policy which instantaneously brings a patient to an anticoagulated state but which requires subsequent infusion rates about 1.5 times greater than I* = .968. By solving the state equations 4, 5, and 11 subject to equations 13 and 14 it is pos­ sible to calculate the infusion policy necessary to keep a patient anticoagulated after the initial time optimal period. Indeed it is this latter infusion period where Ε = 1 and equation 14 holds which determines the 5.25 minute duration of the time optimal per­ iod. As I* becomes larger the duration of the time optimal period becomes shorter and it eventually becomes 0 (solution given by McAvoy {$)). At 5.25 minutes I jumps to I* and then begins to taper in an exponential-like manner as shown in Figure 3. As can be seen, I tapers over a period of 500 minutes and eventually reaches a steady state value of .257 units/hr ml plasma. Also shown in Figure 3 is the response of Ε to the optimal infusion policy. It is interesting to compare the policy shown in Figure 3 to current clinical practice since reasonable parameter values were used in the heparin model. Very often it is found that a patient's need for heparin is greater during the early stages of therapy than later on. In the medical literature this has been referred to as initial resistance to heparin therapy. Generally heparin infusion rates are lowered over a period of several days. The in­ fusion policy shown in Figure 3 does exhibit a tapering transient but the time involved, 500 minutes, is shorter than that observed clinically. This discrepancy may indicate the necessity to refine the model possibly by having the heparin response tied to the slow­ er response of antithrombin III. f

Downloaded by EAST CAROLINA UNIV on January 2, 2018 | http://pubs.acs.org Publication Date: May 30, 1980 | doi: 10.1021/bk-1980-0124.ch024

APPLICATIONS

f

Clinical Use of the Model In presenting the pharmacokinetic model for heparin some parameters were fitted to data on healthy volunteers while other parameters were arbitrarily chosen. To improve heparin therapy the pharmacokinetic modeling approach discussed here must be used in a prospective, predictive manner. By taking additional blood samples early on in heparin therapy data for each patient can be generated. One would measure such variables as heparin level, AT-III level, PTT, etc. Then, by using the model presented here or a more accurate model, parameters can be fitted (7_) and the

Squires and Reklaitis; Computer Applications to Chemical Engineering ACS Symposium Series; American Chemical Society: Washington, DC, 1980.

Squires and Reklaitis; Computer Applications to Chemical Engineering ACS Symposium Series; American Chemical Society: Washington, DC, 1980.

Figure 3.

Optimum dosing policy and clotting time response for heparin model. ( dosing policy; ( ) relative extension of clotting time

) optimum

Downloaded by EAST CAROLINA UNIV on January 2, 2018 | http://pubs.acs.org Publication Date: May 30, 1980 | doi: 10.1021/bk-1980-0124.ch024

Downloaded by EAST CAROLINA UNIV on January 2, 2018 | http://pubs.acs.org Publication Date: May 30, 1980 | doi: 10.1021/bk-1980-0124.ch024

434

COMPUTER

APPLICATIONS

T O CHEMICAL

ENGINEERING

model tailor-made to each individual patient, The present model contains 8 parameters. For prospective clinical useitwould be impractical to try to adaptively f i t this many parameters. No doubt it would be necessary to set many of the parametersinthe model to average constant values. With a small number of remaining parameters it would be possible to adaptively fit these parameters to a patient's data. In this way a heparin model and an optimum dosing policy can be tailor-made for each patient. As more data become available they can be incorporated to further refine the model and in this fashion allow for adaptive control of dosage regimens. An alternative to simply fitting parameters in the model is to attempt to correlate them with easily measured clinical variables. Such variables include: age, height, sex, hematocrit, serum creatinine, etc. Jelliffe and his co-workers have successfully used such a correlation approach to improve the therapeutic results achieved with a number of drugs. For example, they reported (12) on two groups of patients treatedatCarmel Hospital with the anti-arrythmic drug lidocaine. A group of 68 patients constituted the standard therapy group. Of these patients 8 developed ventricular fibrillation, 5 had a toxic reaction and 33 required additional lidocaine. Of 71 patients who received the drug based on a computerized correlation approach, 2 developed ventricular fibrillation, 1 had a toxic reaction and only 2 required additional lidocaine. Similar success has also been reported for digitalis drugs (13). While the above modeling approach looks promising,itis necessary to put it into perspective. The approach is not aimed at having a computer replace a physicianindeciding dosing policy. Rather the approach is proposed as an additional, useful tool which physicians can use to improve therapy. The final decision on any dosing policy must of necessity rest with the attending physician. One way of viewing a physician's role in giving drugs is that heisa control engineer dealing with dynamic processes. The above approach simply formalizes this concept. Conclusions This paper has treated modeling, parameter identification and control of the administration of heparin. A model has been developed using a single CSTR with a volume equal to a patient's blood volume. By using available laboratory data, parameters in the drug model can be fitted and the model tailor-made to each individual patient. Once an accurate model is available it can be used to predict optimum dosing policies. The methods presented have the potential of greatly increasing the efficacy of heparin administration.

Squires and Reklaitis; Computer Applications to Chemical Engineering ACS Symposium Series; American Chemical Society: Washington, DC, 1980.

24. McAVOY Anticoagulant Drug Therapy

435

Abstract A pharmacokinetic model which describes a patient's response to the widely used anticoagulant drug, heparin, is discussed. The model contains several unknown parameters. A patient's own past data can be used to identify these parameters and thus tailor-make a model for each individual under treatment. Once such a model is available it can be used to calculate optimum dosing policies to achieve a specified therapeutic goal.

Downloaded by EAST CAROLINA UNIV on January 2, 2018 | http://pubs.acs.org Publication Date: May 30, 1980 | doi: 10.1021/bk-1980-0124.ch024

Literature Cited 1.

Wessler, S. "Current Dilemmas in the Clinical Use of Heparin," Fed. Proc., 1977, 36, 66-69. 2. Wessler, S. "A Heparin Symposium," Current Ther. Res., 1975, 18, 3-5. 3. Rosenberg, R. "Actions and Interactions of Antithrombin and Heparin," N. Engl. J. Med., 1975, 292, 146-151. 4. Rosenberg, R. "Chemistry of the Hemostatic Mechanism and Its Relationship to the Action of Heparin," Fed. Proc., 1977, 36, 10-19. 5. McAvoy, T. J. "Pharmacokinetic Modeling of Heparin and Its Clinical Implications," J. Pharm. Biopharm., 1979, 7, 331-354. 6. Perrier, D.; Ashley, J.; Levy, G. "Effect of Production Inhibition on Kinetics of Drug Elimination," J. Pharm. Biopharm., 1973, 1, 231-242 7. D'Argenio, D. Z.; Schumitzky, A. "A Program for Simulation and Parameter Estimation in Pharmacokinetic Systems," Comput. Programs in Biomed., 1979, 9, 115-134. 8. Hirsh, J.; VanAken, W. G.; Gallus, S. G.; Dollery, C. T.; Cade, J. F.; Yung, W. L. "Heparin Kinetics in Venous Thrombosis and Pulmonary Embolism," Circulation, 1976,53,691-695, 9. Marciniak, E. "Binding of Heparin in Vitro and In Vivo to Plasma Proteins," J. Lab. Clin. Med., 1974, 84, 344-356. 10. Marciniak, E.; Gockerman, J. "Heparin-Induced Decrease in Circulating Antithrombin III," Lancet, 1977, II, 581-584. 11. Denn, M. "Optimization By Variational Methods"; McGraw-Hill: New York, 1969, pp. 189-191. 12. Jelliffe, R.; Rodman, J.; Kolb, E. "Clinical Studies with Computer-Assisted Lidocaine Infusion Regimens," Circulation, 1976, 54, II-211. 13. Jelliffe, R.; Buell, J.; Kabala, R. "Reduction of Digitalis Toxicity by Computer-Assisted Glycoside Dosage Regimens," Ann. Int. Med., 1972, 77, 891-906. RECEIVED November 5,

1979.

Squires and Reklaitis; Computer Applications to Chemical Engineering ACS Symposium Series; American Chemical Society: Washington, DC, 1980.