Although the therapeutic effects of Vagus Nerve Stimulation (VNS) have been recognized in pre-clinical and pilot clinical studies, the effect of different stimulation configurations on the cardiovascular response is still an open question, especially in the case of VNS delivered synchronously with cardiac activity. In this paper, we propose a formal mathematical methodology to analyze the acute cardiac response to different VNS configurations, jointly considering the chronotropic, dromotropic and inotropic cardiac effects. A latin hypercube sampling method was chosen to design a uniform experimental plan, composed of 75 different VNS configurations, with different values for the main parameters (current amplitude, number of delivered pulses, pulse width, interpulse period and the delay between the detected cardiac event and VNS onset). These VNS configurations were applied to 6 healthy, anesthetized sheep, while acquiring the associated cardiovascular response. Unobserved VNS configurations were estimated using a Gaussian process regression (GPR) model. In order to quantitatively analyze the effect of each parameter and their combinations on the cardiac response, the Sobol sensitivity method was applied to the obtained GPR model and inter-individual sensitivity markers were estimated using a bootstrap approach. Results highlight the dominant effect of pulse current, pulse width and number of pulses, which explain respectively 49.4%, 19.7% and 6.0% of the mean global cardiovascular variability provoked by VNS. More interestingly, results also quantify the effect of the interactions between VNS parameters. In particular, the interactions between current and pulse width provoke higher cardiac effects than the changes on the number of pulses alone (between 6 and 25% of the variability). Although the sensitivity of individual VNS parameters seems similar for chronotropic, dromotropic and inotropic responses, the interacting effects of VNS parameters provoke significantly different cardiac responses, showing the feasibility of a parameter-based functional selectivity. These results are of primary importance for the optimal, subject-specific definition of VNS parameters for a given therapy and may lead to new closed-loop methods allowing for the optimal adaptation of VNS therapy through time.