Spinal cord stimulation is a promising therapy for the recovery after spinal cord injury and control of movement. However, the complex interaction between electrodes and neural reaction impedes more effective neuromodulation. Empirical knowledge and manual adjustment are playing dominant role for the tuning of neuromodulation parameters. The growth of electrode array’s complexity also brings new challenges to find optimal stimulating configurations. Developing a more efficient way to predict the stimulation results and explore the optimal therapy is necessary for the tuning of sensitive neural circuits in spinal cord. We developed a hybrid computational model of human spinal cord and stimulating electrode array to simulate the effect of spinal cord stimulation. The high-precision finite element model of the spinal cord we built matches the anatomical size and segmental innervation of human spinal cord. Biophysical neuron models were embedded into the finite element model. Based on the results of simulation, we analyzed the selectivity of neural stimulation for different muscles. The influence of arrangement and size of the neural electrodes was calculated to provide information for the design of electrode arrays. The polarity configuration and implanting position of the multiple-contact electrode were also considered. We developed Bayesian optimization method to explore the optimal configuration of polarity based on our hybrid model. Our method can efficiently optimize stimulating parameters from a large input space. These results provide effective guidance for electrode design, surgical implantation and neuromodulation therapies for spinal cord stimulation. The optimization result for configuration of polarity shows the potential of using Bayesian optimization in clinical practice. Analysis based on our simulation also contributes to the quantitative understanding of the mechanism of spinal cord stimulation.