Burst suppression is a phenomenon in which the electroencephalogram (EEG) of a pharmacologically sedated patient alternates between higher frequency and amplitude bursting to lower frequency and amplitude suppression. The level of sedation can be quantified by the burst suppression ratio (BSR), which is defined as the amount of time that an EEG is suppressed over the amount of time measured. Maintaining a stable BSR in patients is an important clinical problem, which has led to an interest in the development of BSR-based closed-loop pharmacological control systems. Methods to address the problem typically involve pharmacokinetic (PK) modeling that describes the dynamics of drug infusion in the body as well as signal processing methods for estimating burst suppression from EEG measurements. In this regard, simulations, physiological modeling, and control design can play a key role in producing a solution. This paper seeks to add to prior work by incorporating a Schnider PK model with the Wilson–Cowan neural mass model to use as a basis for robust control design of biophysical burst suppression dynamics. We add to this framework actuator modeling, real-time burst suppression estimation, and uncertainty modeling in both the patient's physical characteristics and neurological phenomena to form a closed-loop system. Using the Ziegler–Nichols tuning method for proportional-integral-derivative (PID) control, we were able to control this system at multiple BSR set points over a simulation time of 18 h in both nominal, patient varied with noise added and with reduced performance due to BSR bounding when including patient variation and noise as well as neurological uncertainty.