Propofol and fentanyl are commonly used agents for the induction of anesthesia, and are often associated with hemodynamic disturbances. Understanding pharmacodynamic impacts is vital for parasympathetic and sympathetic tones during the anesthesia induction period. Inspired by the thermodynamic interaction between drug concentrations and effects, we established a machine-learning-based response surface model (MLRSM) to address this predicament. Then, we investigated and modeled the biomedical phenomena in the autonomic nervous system. Our study prospectively enrolled 60 patients, and the participants were assigned to two groups randomly and equally. Group 1 received propofol first, followed by fentanyl, and the drug sequence followed an inverse procedure in Group 2. Then, we extracted and analyzed the spectrograms of electrocardiography (ECG) and pulse photoplethysmography (PPG) signals after induction of propofol and fentanyl. Eventually, we utilized the proposed MLRSM to evaluate the relationship between anesthetics and the integrity/balance of sympathetic and parasympathetic activity by employing the power of high-frequency (HF) and low-frequency (LF) bands and PPG amplitude (PPGA). It is worth emphasizing that the proposed MLRSM exhibits a similar mathematical form to the conventional Greco model, but with better computational performance. Furthermore, the MLRSM has a theoretical foundation and flexibility for arbitrary numbers of drug combinations. The modeling results are consistent with the previous literature. We employed the bootstrap algorithm to inspect the results’ consistency and measure the various statistical fluctuations. Then, the comparison between the modeling and the bootstrapping results was used to validate the statistical stability and the feasibility of the proposed MLRSM.