Home 
Controlled infusion of intravenous cardiac drugs using global optimization Correspondence Address: OBJECTIVES: The objective of the study is to develop an automatic drug infusion control system during cardiovascular surgery. MATERIALS AND METHODS: Based on the clinical drug dosage analysis, the modeling of cardiovascular system with baroreceptor model is mathematically modeled using compartmental approach, considering the relationship between the volume and flow rate of blood during each heartbeat. This model is then combined with drug modeling of noradrenaline and nitroglycerine by deriving the volume and drug mass concentration equations, based on pharmacokinetics and pharmacodynamics of the drugs. The closedloop patient models are derived from the openloop data obtained from the physiologydrug model with covariate as age. The proportionalintegral controller is designed based on optimal values obtained from bacterial foragingoriented particle swarm optimization algorithm. The controllers are implemented individually for each control variable such as aortic pressure and cardiac output (CO), irrespective of varying weights based on the relative gain array analysis which depicts the maximum influence of cardiac drugs on control variables. RESULTS: The physiologydrug model output responses are simulated using MATLAB. The controlled responses of aortic pressure and CO with infusion rate of cardiac drugs are obtained. The robustness of the controller is checked by introducing variations in cardiovascular model parameters. The efficiency of the controller during normal and abnormal conditions is compared using time domain analysis. CONCLUSIONS: The controller design was efficient and can be further improved by designing switchingbased controllers.
Introduction The operating environment is a highly vulnerable surgical space where each and every patient is going through medical jeopardy. The patients undergoing the surgical regimen are at the advanced stage of their respective ailments. The initial step of the surgical procedure is the drug interventions that are substantially invasive with dynamic and uncertain physiological responses. The intricacy level with respect to cognitive and taskoriented demands, in the surgical suite, results in influential and unrelenting ambiance, which may aggravate the effects of minute human errors and negligence.[1] Surgical process and safety are impacted by various factors such as clinical competence, attributes of the surgical surroundings, surgical workflow, teamwork, and organizational structures. These factors influence each other and any change or malfunction affects the entire system. Coronary artery bypass grafting (CABG) surgery is used to treat and improve the blood flow into the heart. This coronary heart disease is caused by the accumulation of substance called plaque inside the arteries. This plaque may be made of cholesterol, fat found in the blood.[2] The buildup plaque in the walls of arteries causes narrowing which reduces the blood flow into the heart, resulting in heart attack, chest pain, or discomfort. During the surgery, the plagued artery is bypassed or grafted to the healthy artery and this provides a new passage for oxygenrich blood to flow through the heart muscles. The realtime physiological variables and their respective drug infusion rate were obtained where the administration of drug was carried out by an anesthetist manually using syringe pump. The physiological variables obtained during CABG were aortic pressure and cardiac output (CO) with administration of cardiac drugs, such as noradrenaline which acted as vasoconstrictor and nitroglycerine which acted as vasodilator.[3] The maximum allowable limit for noradrenaline is 0.01–3 μg/kg and for nitroglycerine is 0.01–1.5 μg/kg in clinical practice. Materials and Methods Mathematical modeling Cardiovascular modeling The physiological variable such as aortic pressure from which the mean arterial pressure (MAP) is manipulated and CO is controlled using the controller designed by regulating the infusion rate of cardiac drugs.[4] The controller design is based on the mathematical model of the process involved. The mathematical model consists of cardiovascular system with baroreflex mechanism combined with drug modeling based on pharmacokinetics and pharmacodynamics.[5] Since these drugs are used to regulate the hemodynamics of the heart, the cardiovascular model is initially derived by considering four compartments, namely left heart, systemic circulation, right heart, and pulmonary circulation.[6] The circulatory compartments are partitioned into three arterial and two ventricular sections. When the left ventricular pressure exceeds the systemic aortic pressure, the aortic valve opens and the blood enters into systemic circulation through peripheral systemic resistance.[7] [INLINE:1] Q1v = 0, when aortic value is closed [INLINE:2] [INLINE:3] [INLINE:4] Where Q1v is left ventricular flow rate Q1v is inertial vessel properties of left ventricle P1v is left ventricular pressure Pas is systemic aortic pressure V1v is left ventricular volume E1v is elastance of left ventricle Vd1v is left ventricular volume at zero pressure Rsys is systemic resistance Pa1 is pressure at systemic arterial section 1. Systemic circulatory arterial section 1 [INLINE:5] [INLINE:6] [INLINE:7] Systemic circulatory arterial section 2 [INLINE:8] [INLINE:9] Systemic circulatory arterial section 3 [INLINE:10] [INLINE:11] Where Qa1 is the flow rate of arterial section 1 Ra1, Ra2, Ra3 are the viscosity of arterial sections 1, 2, and 3 Ca1, Ca2, Ca3 are the elastic properties of arterial sections 1, 2, and 3 La1 is the inertial vessel properties of arterial section 1 Pa2, P a3 are the pressure at arterial sections 2 and 3 Va1, Va2, Va3, are the volume of arterial sections 1, 2, and 3 Vuna1 is the volume of arterial section 1 at zero pressure Vuna3 is the volume of arterial section 3 at zero pressure. Systemic circulatory venous section 1 [INLINE:12] [INLINE:13] Systemic circulatory venous section 2 [INLINE:14] [INLINE:15] [INLINE:16] Where Qv2 is the flow rate of venous section 1 Rv1, Rv2 are the viscosity of venous sections 1 and 2 Cv1, Cv2 are the elastic properties of venous section 1 and 2 Lv2 is the inertial vessel properties of venous section 2 Pv1, Pv2 are the pressure at venous sections 1 and 2 Vv1, Vv2 are the volume of venous sections 1 and 2 Vunv1 is the volume of venous section 1 at zero pressure Vunv2 is the volume of venous section 2 at zero pressure. The veins return the blood to right atrium. When the right atrial pressure exceeds right ventricular pressure, the tricuspid valve opens and the blood enters into right ventricle. When the right ventricular pressure exceeds the pulmonary aortic pressure, the blood enters into pulmonary circulation through pulmonary valve. [INLINE:17] Qrv = 0, when pulmonary valve is closed [INLINE:18] [INLINE:19] [INLINE:20] Where Qrv is right ventricular flow rate Lrv is inertial vessel properties of right ventricle Prv is right ventricular pressure Pap is aortic pulmonary pressure Vrv is right ventricular volume Erv is elastance of right ventricle Vdrv is right ventricular volume at zero pressure Rpul is pulmonary resistance Pp1 is pressure at pulmonary arterial section 1. Pulmonary circulatory arterial section 1 [INLINE:21] [INLINE:22] [INLINE:23] Pulmonary circulatory arterial section 2 [INLINE:24] [INLINE:25] Pulmonary circulatory arterial section 3 [INLINE:26] [INLINE:27] Where Qp1 is the flow rate of pulmonary section 1 Rp1, Rp2, Rp3 are the viscosity of pulmonary section 1, 2, and 3 Cp1, Cp2, Cp3 are the elastic properties of pulmonary sections 1, 2, and 3 Lp1 is the inertial vessel properties of pulmonary section 1 Pp2, Pp3 are the pressure at pulmonary sections 2 and 3 Vp1, Vp2, Vp3 are the volume of pulmonary sections 1, 2, and 3 Vunp1 is the volume of pulmonary section 1 at zero pressure Vunp3 is the volume of pulmonary section 3 at zero pressure. Pulmonary circulatory venous section 1 [INLINE:28] [INLINE:29] Pulmonary circulatory venous section 2 [INLINE:30] [INLINE:31] [INLINE:32] Where Ql2 is the flow rate of venous section 1 Rl1, Rl2 are the viscosity of venous sections 1 and 2 Cl1, Cl2 are the elastic properties of venous sections 1 and 2 Ll2 is the inertial vessel properties of venous section 2 Pl1, Pl2 are the pressure at venous sections 1 and 2 Vl1, Vl2 are the volume of venous sections 1 and 2 Vunl1 is the volume of venous section 1 at zero pressure Vunl2 is the volume of venous section 2 at zero pressure. The veins return the blood to the left atrium. When the left atrial pressure exceeds the left ventricular pressure, the mitral valve opens and the blood enters into left ventricle. In this way, the cardiovascular system is modeled which represents the continuous closed loop system with aortic pressure and CO passing into baroreceptor model. [INLINE:33] Qla = 0 when mitral valve is closed [INLINE:34] [INLINE:35] Where Qla is flow rate of left atria Lla is inertial vessel properties of left atria Pla is left atrial pressure Plv is left ventricular pressure Rla is resistance of left atria Vla is volume of left atria Ql2 is flow rate of pulmonary venous section2 Ela is elastance of left atria Vdla is volume of left atria at zero pressure. Baroreceptor modeling Baroreceptor model acts as the feedback system which helps in shortterm regulation of blood pressure. Baroreceptor is stretch receptor located in the walls of the blood vessels, the most accessible of which are located in the carotid sinus and in the aortic arch.[8] Carotid sinus baroreceptor is located in the distinctive part of the two common carotid sinus arteries. Aortic arch baroreceptor, on the other hand, is located in the walls of the aortic arch. The aortic arch and the carotid sinus receptors are believed to be functionally same, expect that the aortic arch receptors operate at a higher pressure level.[9] The baroreceptor modeled here is concentrated on the carotid sinus which has nerve ending and responds to deformation in the walls of the blood vessels. The nerve activity evolves from two components, a pressuremechanical and mechanicalelectrical component. The nerve activity from the carotid sinus receptor is called firing rates. The baroreceptor acts immediately to deformations in the blood vessel walls by a pressuremechanical mechanism. The mechanicalelectrical component generates the baroreceptor firing rate by the mechanicalelectrical mechanism inbuilt in the receptor.[10] The cardioinhibitory center and vasomotor center of the central nervous system generates sympathetic nerve activity. The efferent pathway transmits the nerve activity to the cardiovascular model. The input function to the baroreceptor model is expressed as: [INLINE:36] Where MAPnom is the nominal value of MAP and α is the constant. The parameters that are influenced by autonomic reflexes are modeled in such a way that it provides continuous control action to regulate the MAP as shown. [INLINE:37] [INLINE:38] [INLINE:39] [INLINE:40] Where [INSIDE:1] represents the output from baroreceptor model for that instant. The nominal values such as [INSIDE:2]. The aortic pressure and CO from cardiovascular and baroreceptor model are shown in [Figure 1].{Figure 1} Pharmacokineticpharmacodynamic modeling The pharmacokinetics and pharmacodynamics of the cardiac drugs such as noradrenaline and nitroglycerine are studied by deriving their volume and drug mass concentration equation which is intervened into the blood and flows into cardiovascular system through circulation process.[11] The flow relationship in four compartments are modeled as: [INLINE:41] Where x represents the compartment being considered, [INSIDE:3] is the input flow from previous compartment and [INSIDE:4] is the output flow from the next compartment. The amount of drug entering into each compartment is modeled as [INLINE:42] Where Concdg= M/V, M is the mass of drug in compartment and V is the volume of the compartment and τ1/2 is the halflife of the drug in the compartment. Pharmacodynamics is completely based on the circulatory parameters with the effect site concentration of cardiac drugs [INLINE:43] Where Eff is the measure of drug effect on the affected parameter in the compartment where the effect is considered. f1 is determined from f1= f2/p. The termpis expressed as: [INLINE:44] Here the total blood volume is assumed to be 85 ml/kg of body weight. f1 and f2 are drug constants, Effmax is the maximum drug effect, b50 represents 50% of drug effect on corresponding infusion, and L is the power of maximum concentration in drug effect mass equation. The effect of noradrenaline and nitroglycerine on control variables are shown in [Figure 2] and [Figure 3].{Figure 2}{Figure 3} Based on the physiologydrug model derived, the open loop data were obtained for patient model with weight of 65 kg by varying the infusion rate within the allowable limit. The MAP is calculated from aortic pressure using the formulae: MAP = ⅓systolic pressure + ⅔diastolic pressure The transfer function of the patient model was determined using autoregressive model with exogenous input (ARX) technique. The patient model is a 2 × 2 system represented as: [INLINE:45] Since the system derived is a multiinputmultioutput system, the maximum influence of cardiac drugs on control variables is to be determined. This is done using relative gain array (RGA) analysis. RGA uses the steady gain of the process to calculate the drug interactions using the expression [INLINE:46] k11, k12, k21, and k22 are the steady state gain of the 2 × 2 process. The RGA value obtained for the patient model is 1.000000049 with noradrenaline having maximum influence on CO and nitroglycerine having maximum influence on MAP. Bacterial foragingoriented particle swarm optimization Optimization techniques are involved in process with large complexity and nonlinearity to estimate the optimal values to control the process.[12],[13],[14],[15] In this study, the proportionalintegral (PI) controller tuning parameters are determined using bacterial foragingoriented particle swarm optimization (BFOAPSO). The combination of bacterial foraging and particle swarm techniques are implemented here.[16] The BFOA and PSO algorithms are combined to utilize the ability of PSO algorithm for social information tradeoff and ability of BFOA to determine new solution using elimination and dispersal stage. The bacterial foraging algorithm consists of four stages – chemotaxis, swarming, reproduction, and eliminationdispersal.[17] The MATLAB coding for BFOAPSO algorithm is as follows: %% Tuning of PI controller using bacterial foraging oriented particle swarm optimization %Initialization clear all clc d = 2; % dimension of search space Tb = 10; % The number of bacteria Nch = 5; % Number of chemotactic steps Nsl = 4; % Limits the length of a swim Nrep = 2; % The number of reproduction steps Ned = 2; % The number of eliminationdispersal events Sr = s/2; % The number of bacteria reproductions per generation Ped = 0.25; % The probability that each bacteria will be eliminated/dispersed x(:,1)=0.001*ones (Tb, 1); % the run length for h = 1:Tb Delta(:, h)=(2*round (rand (p, 1))1).*rand (p, 1); end for k = 1:Tb% the initial positions P (1,:,1,1,1)=50*rand (Tb, 1)'; P (2,:,1,1,1)=0.2*rand (Tb, 1)'; %P (3,:,1,1,1)=0.2*rand (Tb, 1)'; end c1 = 1.2; % PSO parameter C1 1.2 c2 = 0.5; % PSO parameter C2.5 R1 = rand(p, s); % PSO parameter R2 = rand(p, s); % PSO parameter Plocal_best_position = 0*ones (p, Tb, Nch); % PSO Pglobal_best_position = 0*ones (p, Tb, Nch); % PSO velocity =0.3*randn (p, Tb); % PSO current_position= 0*ones (p, Tb, Nch); % PSO %Main loop %Elimination and dispersal loop for v=1:Ned %Reproduction loop for m=1:Nrep %Swim/tumble (chemotaxis) loop for l=1:Nch for h=1:Tb J(h, l, m, v)=tracklsq (P(:, h, l, m, v)); % Tumble Jlast= J (h, l, m, v); Jlocal (h, l)=Jlast; P(:, h, l+1, m, v)=P(:, h, l, m, v)+x (h, m)*Delta(:, h); % This adds a unit vector in the random direction % Swim (for bacteria that seem to be headed in the right direction) J (h, l+1, m, v)=tracklsq (P(:, h, l+1, m, v)); k=0; % Initialize counter for swim length while ksl k = k+1; if J (h, l+1, m, v)last Jlast= J (h, l+1, m, v); P(:, h, l+1, m, v)=P(:, h, l+1, m, v)+x (h, m)*Delta (:, h); J (h, l+1, m, v)=tracklsq (P(:, h, l+1, m, v)); current_position(:, h, l+1)= P(:, h, l+1, m, v); % PSO Jlocal(h, l+1) = J (h, l+1, m, v); % PSO else Jlocal(h, l+1) = J (h, l+1, m, v); current_position(:, h, l+1)= P(:, h, l+1, m, v); k = Nsl; end end sprintf('The value of interation h %3.0f, l = %3.0f, m= %3.0f, v= %3.0f', h, l, m, v); end % Go to next bacterium %For each chemotactic evaluate the local and global best position %Local best position over all chemotactic that each bacteria move through it [Jmin_for_each_chemotactic, index]=min (Jlocal,[],2); for k = 1:Tb Plocal_best_position(:, k, l) = current_position(:, k, index (k,:)); end %Global best position over all chemotactic and for each bacteria [Y, I]=min (Jmin_for_each_chemotactic); global_best_position = current_position(:, I, index (I,:)); for k = 1:Tb Pglobal_best_position(:, k, l)=global_best_position; end %Calculate the new direction for each bacteria velocity =0.9* velocity + c1*(R1.*(Plocal_best_position(:,:, l)current_position(:,:l+1))) + c2*(R2.*(Pglobal_best_position(:,:, l)current_position(:,:, l+1))); Delta = velocity; end % Go to the next chemotactic %Reprodution Jhealthy= sum (J(:,:, m, v),2); % Set the health of each of the S bacteria [Jhealthy, sortind]=sort (Jhealthy); % Sorts the nutrient concentration in order of ascending P(:,:,1, m+1, v)=P(:, sortind, Nch+1, m, v); x(:, m+1)=x (sortind, m); % And keeps the chemotaxis parameters with each bacterium at the next generation %Split the bacteria (reproduction) for h=1:Sr P(:, h+Sr, 1, m+1, v)=P(:, h, l, m+1, v); % The least fit do not reproduce, the most fit ones split into two identical copies x (h+Sr, m+1)=x (h, m+1); end end % Go to next reproduction %Elimination and dispersal for k=1:Tb if Ped>rand % % Generate random number P (1,:,1,1,1)=0.2*rand (s, 1)'; P (2,:,1,1,1)=0.2*rand (s, 1)'; %P (3,:,1,1,1)=0.2*rand (s, 1)'; else P(:, k, 1,1, v+1)=P(:, k, 1, Nrep+1, v); % Bacteria that are not dispersed end end end % Go to next elimination and dispersal %Report reproduction = J(:,1:Nch, Nrep, Ned); [jlastreproduction, O] = min (reproduction,[],2); % min cost function for each bacterial [Y, I]=min (jlastreproduction); pbest=P(:, I, O (I,:), m, v); Kp=pbest (1,:) Ki=pbest (2,:) Results The clinical data observed during CABG surgery are obtained as follows. The variations in MAP and CO with changes in noradrenaline and nitroglycerine infusions are shown in [Figure 4], [Figure 5], [Figure 6]. The maximum allowable limit for noradrenaline infusion is 0.01–3 μg/kg and nitroglycerine infusion is 0.01–1.5 μg/kg. The optimal values for PI controller were obtained and applied in patient model to simulate the controlled responses. [Figure 7] shows the controlled MAP and CO with infusion of noradrenaline and nitroglycerine. The infusion rates of noradrenaline and nitroglycerine are shown in [Figure 8]. The controller sensitivity was tested by introducing a disturbance at 500 s for MAP and CO. Since the MAP is influenced by nitroglycerine infusion, MAP is in decreasing manner from 150 mmHg to desired MAP of 93 mmHg and CO is influenced by infusion of noradrenaline hence the blood flow is in increasing manner from 3300 ml to 5000 ml. The controlled infusion rate of nitroglycerine is 0.02 μg/kg and that of noradrenaline is 0.1 μg/kg.{Figure 4}{Figure 5}{Figure 6}{Figure 7}{Figure 8} Discussion The robustness of the controller is verified by considering an abnormality in the cardiovascularbaroreceptor model. The abnormality is incorporated by varying the systemic resistance from 0.0334 mmHg*s/ml to 0.1 mmHg*s/ml and the patient model obtained was [INLINE:47] [INLINE:48] The same PI controller is used in this patient model and the robustness of the controller based on optimization algorithm is shown in [Figure 9] and [Figure 10]. The controlled infusion rate of nitroglycerine is 0.02 μg/kg and that of noradrenaline is 0.2 μg/kg as shown in [Figure 10]. The efficiency and performance of the controller is determined using time domain analysis as tabulated in [Table 1].{Figure 9}{Figure 10}{Table 1} Conclusions Infusion of multiple drugs creates a significant amount of complexity in the modeling process due to interaction of drugs. The future work focuses on study about interaction of drugs from clinical data and developing switchingbased controllers for the efficient drug administration. Furthermore, the developed model and their respective controller which provides control action on cardiac drugs and anesthetic agent should be evaluated in the clinical setting to guarantee the total reliability of the controllers developed. The clinical data collected from this clinical evaluation can also be utilized in finetuning of the controllers. Financial support and sponsorship Nil. Conflicts of interest There are no conflicts of interest. References


