Background: A spinal motoneuron contacts a bunch of muscle fibers forming a motor unit that underlies all mammalian movements. The essential role of the motor unit is the transduction of synaptic inputs from descending and reflex pathways into muscle force. Since the input–output properties of both motoneurons and muscle fibers are non-linear, it has been difficult to make predictions on how changes in synaptic inputs to motoneuron, cellular properties of the motoneuron and muscle fibers and muscle length may affect motor output [1]. Methods: To tackle this fundamental issue in the field of motor neuroscience, we developed a physiologically plausible but computationally efficient model of the motor unit and a software package that allows for virtual experiments on the input–output properties of the motor unit over a full range of physiological inputs and biophysical parameters. Results: The computational model of motor unit was first built in this study coupling the motoneuron model and the muscle unit model with a simplified axon model. The motoneuron model was developed using the recently reported two-compartment modeling approach [2]. The key feature of the new reduced motoneuron model is that all cable parameters of the reduced model are analytically determined based on the system properties such as input resistance, membrane time constant and electrical coupling properties between the soma and the dendrites, which are all empirically measurable from real motoneurons. For the muscle unit, the recently developed muscle modeling approach was employed that consists of three sub-modules representing [3]: (1) the transformation of the spike signals from motoneurons into the dynamics of calcium concentration in the sarcoplasm, (2) the conversion of the calcium concentration to the muscle activation level, and (3) the transformation of the muscle activation level into the muscle force using Hill-type muscle mechanics. The new muscle model was constructed in this study to reflect all experimentally identified dependencies of muscle activation dynamics on muscle length and movement over a full range of stimulation frequencies in cat soleus muscles. Then, to enhance the usability and extendibility the software package for simulating and analyzing the developed motor unit model was designed and implemented based on the object-oriented programing paradigm and open source Python language along with graphic user interfaces (GUI). The software package developed in this study provides a GUI-based simulation environment in which a single motoneuron, muscle unit, and motor unit can be individually simulated and analyzed in a wide range of experimental conditions reflecting biological realisms. Conclusions: Our model of the motor unit and user-friendly simulation software may provide not only a computational framework to gain systemic insights into motor control by the central nervous system in a cellular perspective but also a basis on which to build biologically realistic large-scale neuro-musculo-skeletal models.