

RESEARCH ARTICLE 



Year : 2019  Volume
: 51
 Issue : 1  Page : 6171 

Controlled infusion of intravenous cardiac drugs using global optimization
GC Sowparnika^{1}, M Thirumarimurugan^{1}, VM Sivakumar^{1}, N Vinoth^{2}
^{1} Department of Chemical Engineering, Coimbatore Institute of Technology, Coimbatore, Tamil Nadu, India ^{2} Department of Instrumentation Engineering, Madras Institute of Technology, Chennai, Tamil Nadu, India
Date of Submission  04Feb2019 
Date of Acceptance  20Feb2019 
Date of Web Publication  19Mar2019 
Correspondence Address: Ms. G C Sowparnika Department of Chemical Engineering, Coimbatore Institute of Technology, Coimbatore  641 014, Tamil Nadu India
Source of Support: None, Conflict of Interest: None  Check 
DOI: 10.4103/ijp.IJP_612_18
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.
Keywords: Cardiac drugs, infusion rate, modeling, optimization, time domain analysis
How to cite this article: Sowparnika G C, Thirumarimurugan M, Sivakumar V M, Vinoth N. Controlled infusion of intravenous cardiac drugs using global optimization. Indian J Pharmacol 2019;51:6171 
How to cite this URL: Sowparnika G C, Thirumarimurugan M, Sivakumar V M, Vinoth N. Controlled infusion of intravenous cardiac drugs using global optimization. Indian J Pharmacol [serial online] 2019 [cited 2020 May 25];51:6171. Available from: http://www.ijponline.com/text.asp?2019/51/1/61/254589 
» 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]}
Q_{1v} = 0, when aortic value is closed
Where Q_{1v} is left ventricular flow rate
Q_{1v} is inertial vessel properties of left ventricle
P_{1v} is left ventricular pressure
P_{as} is systemic aortic pressure
V_{1v} is left ventricular volume
E_{1v} is elastance of left ventricle
V_{d1v} is left ventricular volume at zero pressure
R_{sys} is systemic resistance
P_{a1} is pressure at systemic arterial section 1.
Systemic circulatory arterial section 1
Systemic circulatory arterial section 2
Systemic circulatory arterial section 3
Where Q_{a1} is the flow rate of arterial section 1
R_{a1}, R_{a2}, R_{a3} are the viscosity of arterial sections 1, 2, and 3
C_{a1}, C_{a2}, C_{a3} are the elastic properties of arterial sections 1, 2, and 3
L_{a1} is the inertial vessel properties of arterial section 1
P_{a2}, P _{a3} are the pressure at arterial sections 2 and 3
V_{a1}, V_{a2}, V_{a3}, are the volume of arterial sections 1, 2, and 3
V_{una1} is the volume of arterial section 1 at zero pressure
V_{una3} is the volume of arterial section 3 at zero pressure.
Systemic circulatory venous section 1
Systemic circulatory venous section 2
Where Q_{v2} is the flow rate of venous section 1
R_{v1}, R_{v2} are the viscosity of venous sections 1 and 2
C_{v1}, C_{v2} are the elastic properties of venous section 1 and 2
L_{v2} is the inertial vessel properties of venous section 2
P_{v1}, P_{v2} are the pressure at venous sections 1 and 2
V_{v1}, V_{v2} are the volume of venous sections 1 and 2
V_{unv1} is the volume of venous section 1 at zero pressure
V_{unv2} 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.
Q_{rv} = 0, when pulmonary valve is closed
Where Q_{rv} is right ventricular flow rate
L_{rv} is inertial vessel properties of right ventricle
P_{rv} is right ventricular pressure
P_{ap} is aortic pulmonary pressure
V_{rv} is right ventricular volume
E_{rv} is elastance of right ventricle
V_{drv} is right ventricular volume at zero pressure
R_{pul} is pulmonary resistance
P_{p1} is pressure at pulmonary arterial section 1.
Pulmonary circulatory arterial section 1
Pulmonary circulatory arterial section 2
Pulmonary circulatory arterial section 3
Where Q_{p1} is the flow rate of pulmonary section 1
R_{p1}, R_{p2}, R_{p3} are the viscosity of pulmonary section 1, 2, and 3
C_{p1}, C_{p2}, C_{p3} are the elastic properties of pulmonary sections 1, 2, and 3
L_{p1} is the inertial vessel properties of pulmonary section 1
P_{p2}, P_{p3} are the pressure at pulmonary sections 2 and 3
V_{p1}, V_{p2}, V_{p3} are the volume of pulmonary sections 1, 2, and 3
V_{unp1} is the volume of pulmonary section 1 at zero pressure
V_{unp3} is the volume of pulmonary section 3 at zero pressure.
Pulmonary circulatory venous section 1
Pulmonary circulatory venous section 2
Where Q_{l2} is the flow rate of venous section 1
R_{l1}, R_{l2} are the viscosity of venous sections 1 and 2
C_{l1}, C_{l2} are the elastic properties of venous sections 1 and 2
L_{l2} is the inertial vessel properties of venous section 2
P_{l1}, P_{l2} are the pressure at venous sections 1 and 2
V_{l1}, V_{l2} are the volume of venous sections 1 and 2
V_{unl1} is the volume of venous section 1 at zero pressure
V_{unl2} 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.
Q_{la} = 0 when mitral valve is closed
Where Q_{la} is flow rate of left atria
L_{la} is inertial vessel properties of left atria
P_{la} is left atrial pressure
P_{lv} is left ventricular pressure
R_{la} is resistance of left atria
V_{la} is volume of left atria
Q_{l2} is flow rate of pulmonary venous section2
E_{la} is elastance of left atria
V_{dla} 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:
Where MAP_{nom} 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.
Where represents the output from baroreceptor model for that instant. The nominal values such as . The aortic pressure and CO from cardiovascular and baroreceptor model are shown in [Figure 1].  Figure 1: Aortic pressure (a) and cardiac output (b) obtained from cardiovascularbaroreceptor model
Click here to view 
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:
Where x represents the compartment being considered, is the input flow from previous compartment and is the output flow from the next compartment. The amount of drug entering into each compartment is modeled as
Where Conc_{dg}= 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
Where Eff is the measure of drug effect on the affected parameter in the compartment where the effect is considered. f1 is determined from f_{1}= f_{2}/p. The termpis expressed as:
Here the total blood volume is assumed to be 85 ml/kg of body weight. f_{1} and f_{2} are drug constants, Eff_{max} 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: Noradrenaline effect on aortic pressure (a) and cardiac output (b)
Click here to view 
 Figure 3: Nitroglycerine effect on aortic pressure (a) and cardiac output (b)
Click here to view 
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:
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
k_{11}, k_{12}, k_{21}, and k_{22} 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
T_{b} = 10; % The number of bacteria
N_{ch} = 5; % Number of chemotactic steps
N_{sl} = 4; % Limits the length of a swim
N_{rep} = 2; % The number of reproduction steps
N_{ed} = 2; % The number of eliminationdispersal events
S_{r} = s/2; % The number of bacteria reproductions per generation
P_{ed} = 0.25; % The probability that each bacteria will be eliminated/dispersed
x(:,1)=0.001*ones (T_{b}, 1); % the run length
for h = 1:T_{b}
Delta(:, h)=(2*round (rand (p, 1))1).*rand (p, 1);
end
for k = 1:T_{b}% the initial positions
P (1,:, 1, 1,1)=50*rand (T_{b}, 1)';
P (2,:, 1, 1,1)=0.2*rand (T_{b}, 1)';
%P (3,:, 1, 1,1)=0.2*rand (T_{b}, 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
P_{local_best_position} = 0*ones (p, T_{b}, N_{ch}); % PSO
P_{global_best_position} = 0*ones (p, T_{b}, N_{ch}); % PSO
velocity =0.3*randn (p, T_{b}); % PSO
current_{_position}= 0*ones (p, T_{b}, N_{ch}); % PSO
%Main loop
%Elimination and dispersal loop
for v=1:N_{ed}
%Reproduction loop
for m=1:N_{rep}
%Swim/tumble (chemotaxis) loop
for l=1:N_{ch}
for h=1:T_{b}
J(h, l, m, v)=tracklsq (P(:, h, l, m, v));
% Tumble
J_{last}= J (h, l, m, v);
J_{loca}l (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
J_{last}= 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
Jl_{ocal}(h, l+1) = J (h, l+1, m, v); % PSO
else
J_{local}(h, l+1) = J (h, l+1, m, v);
current_{_position}(:, h, l+1)= P(:, h, l+1, m, v);
k = N_{sl};
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
[J_{min}__{for_each_chemotactic, index}]=min (J_{local},[],2);
for k = 1:T_{b}
P_{local_best_position}(:, k, l) = current_{_position}(:, k, index (k,:));
end
%Global best position over all chemotactic and for each bacteria
[Y, I]=min (J_{min_for_each_chemotactic});
global_{_best_position} = current_{_position}(:, I, index (I,:));
for k = 1:T_{b}
P_{global_best_position}(:, k, l)=global_best_position;
end
%Calculate the new direction for each bacteria
velocity =0.9* velocity + c1*(R1.*(P_{local_best_position}(:,:, l)current_{_position}(:,:l+1))) + c2*(R2.*(P_{global_best_position}(:,:, l)current_{_position}(:,:, l+1)));
Delta = velocity;
end % Go to the next chemotactic
%Reprodution
J_{healthy}= sum (J(:,:, m, v),2); % Set the health of each of the S bacteria
[J_{healthy}, sortind]=sort (J_{healthy}); % Sorts the nutrient concentration in order of ascending
P(:,:, 1, m+1, v)=P(:, sortind, N_{ch}+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:S_{r}
P(:, h+S_{r}, 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+S_{r}, m+1)=x (h, m+1);
end
end % Go to next reproduction
%Elimination and dispersal
for k=1:T_{b}
if P_{ed}>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, N_{rep}+1, v); % Bacteria that are not dispersed
end
end
end % Go to next elimination and dispersal
%Report
reproduction = J(:,1:N_{ch}, N_{rep}, N_{ed});
[j_{lastreproduction}, O] = min (reproduction,[],2); % min cost function for each bacterial
[Y, I]=min (j_{lastreproduction});
pbest=P(:, I, O (I,:), m, v);
K_{p}=pbest (1,:)
K_{i}=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: Mean arterial pressure observed during coronary artery bypass grafting with (a) noradrenaline and (b) nitroglycerine infusion
Click here to view 
 Figure 5: Cardiac output observed during coronary artery bypass grafting with (a) noradrenaline and (b) nitroglycerine infusion
Click here to view 
 Figure 6: Infusion rate of noradrenaline (a) and nitroglycerine (b) observed during coronary artery bypass grafting surgery
Click here to view 
 Figure 7: Mean arterial pressure (a) and cardiac output (b) obtained using bacterial foragingoriented particle swarm optimizationbased proportionalintegral controller
Click here to view 
 Figure 8: Controlled nitroglycerine (a) and noradrenaline (b) infusion rate
Click here to view 
» 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
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: Mean arterial pressure (a) and cardiac output (b) obtained using bacterial foragingoriented particle swarm optimizationbased proportionalintegral controller for change in systemic resistance
Click here to view 
 Figure 10: Controlled nitroglycerine (a) and noradrenaline (b) infusion obtained using bacterial foragingoriented particle swarm optimizationbased proportionalintegral controller for change in systemic resistance
Click here to view 
» 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   
1.  Anju C, Nafeesa K. Control scheme for arterial blood pressure regulation in hypertensive patients. Int J Adv Inf Sci Technol 2014;30:51720. 
2.  Arpita B, Ashoke S. Online identification and internal model control for regulating hemodynamic variables in congestive heart failure patient. Int J Phar Med Biol Sci 2015;24:859. 
3.  Nirmala SA, Muthu R, Veena Abirami B. Model predictive control of drug infusion system for mean arterial pressure regulation of critical care patients. Res J Appl Sci Eng Technol 2014;7:46015. 
4.  Chen B, Song T, Guo T, Xiang H, Liu Y, Qin Y, et al. Asimplified computer model of cardiovascular system with an arm branch. Biomed Mater Eng 2014;24:255561. 
5.  Marques AG, SogueroRuiz C, Ramos J, MoraJimenez I, GoyaEsteban R, GarciaCarretero R, et al. Modelling cardiovascular condition evolution in hypertensive population using graph signal processing. Comput Cardiol 2017;44:14. 
6.  Sarabadani Tafreshi A, KlamrothMarganska V, Nussbaumer S, Riener R. Realtime closedloop control of human heart rate and blood pressure. IEEE Trans Biomed Eng 2015;62:143442. 
7.  Le TQ, Bukkapatnam ST, Komanduri R. Realtime lumped parameter modeling of cardiovascular dynamics using electrocardiogram signals: Toward virtual cardiovascular instruments. IEEE Trans Biomed Eng 2013;60:235060. 
8.  Hegyi G, Drzewiecki G. Nonlinear Dynamic Model of Baroreceptor Blood Pressure Regulation. Boston, USA: Proceedings of the 40 ^{th} Annual Northeast Bioengineering Conference; 2527 April, 2014. p. 12. 
9.  Ojeda D, Rolle VL, Rossel O, Karam N, Hagege A, Bonnet JL, et al. Analysis of a Baroreflex Model for the Study of Chronotropic Response to Vagal Nerve Stimulation. Montpellier, France: Proceedings of the 7 ^{th} International IEEE/EMBS Conference on Neural Engineering; 2224 April, 2015. p. 5414. 
10.  Simpson DM, Beda A. Individual Difference in Baroreceptor Sensitivity between Increasing and Decreasing Blood Pressure Sequences. Trento, Italy: Proceedings of the 8 ^{th} Conference of European Study Group on Cardiovascular Oscillations; 2528 May, 2014. p. 1612. 
11.  Bighamian R, Soleymani S, Reisner AT, Seri I, Hahn JO. Prediction of hemodynamic response to epinephrine via modelbased system identification. IEEE J Biomed Health Inform 2016;20:41623. 
12.  AbdElazim SM, Ali E. A hybrid particle swarm optimization and bacterial foraging for optimal power system stabilizers design. Int J Electr Power 2013;46:33441. 
13.  Devi S, Geethanjali M. Application of Modified Bacterial Foraging Optimization algorithm for optimal placement and sizing of Distributed Generation. Expert Syst Appl 2014; 41: 277281. 
14.  Kora P, Krishna KS. Hybrid firefly and particle swarm optimization algorithm for the detection of bundle branch block. Indian J Clin Anaesth 2016;2:448. 
15.  Niu B, Wang H, Wang J, Tan L. Multiobjective bacterial foraging optimization. Neurocomputing 2013;116:33645. 
16.  Ozedmir MT, Ozturk D, Eke I, Celik V, Lee KY. Tuning of optimal classical and fractional PID parameters for automatic generation control based on the bacterial swarm optimization. IFAC Pap Online 2015;48:5016. 
17.  Peng CY, Ying L, Gang W, Feng ZY, Qian X, Hao FJ, et al. A novel bacterial foraging optimization algorithm for feature selection. Expert Syst Appl 2017;83:117. 
[Figure 1], [Figure 2], [Figure 3], [Figure 4], [Figure 5], [Figure 6], [Figure 7], [Figure 8], [Figure 9], [Figure 10]
[Table 1]
