The neuregulin/ErbB system is a growth factor/receptor cascade that has been proven to be essential in the development of the heart and the sympathetic nervous system. However, the basis of the specificity of ligand-receptor recognition remains to be elucidated. In this study, the structures of NRG-1beta/ErbB3 and NRG-1beta/ErbB4 complexes were modeled based on the available structures of the homologous proteins. The binding free energies of NRG-1beta to ErbB3 and ErbB4 were calculated using the molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) computational method. In addition, computational alanine-scanning mutagenesis was performed in the binding site of NRG-1beta and the difference in the binding free energies between NRG-1beta mutants and the receptors was calculated. The results specify the contribution of each residue at the interaction interfaces to the binding affinity of NRG-1beta with ErbB3 and ErbB4, identifying several important interaction residue pairs that are in agreement with previously acquired experimental data. This indicates that the presented structural models of NRG-1beta/ErbB3 and NRG-1beta/ErbB4 complexes are reliable and could be used to guide future studies, such as performing desirable mutations on NRG-1beta to increase the binding affinity and selectivity to the receptor and discovering new therapeutic agents for the treatment of heart failure.