Hepatitis B core antigen virus-like particles (HBc-VLPs) are widely used as vaccine vectors due to their good stability and easy modification, and the investigations of VLPs is one of the hot spots in the field of bio-pharmaceutical engineering. However, VLPs may disassemble or aggregate due to their sensitivity to temperature, pH and other factors, which becomes the bottleneck hindering the widely application, and the underlying mechanisms which governs the structure and thus stability of VLPs is still ambiguous. In this work, molecular dynamics simulation was utilized to investigate the stability of the dimer, pentamer and hexamer formed by protein subunits in HBc-VLP. Instead of using empirical values in previous studies, the parameters of protein dielectric constant in aqueous solution were obtained by molecular dynamics simulations, and the results suggested that both the aqueous solvent and the arrangement of protein subunits in the complex could significantly change the dielectric constant, which further affected the binding free energy. Furthermore, with the dielectric constant of protein subunits, the binding free energy between protein subunits were calculated by the molecular mechanics-Poisson Boltzmann solvent accessible surface area (MM-PBSA) method. Finally, according to the calculation results, it was speculated that the stability of the hexamer was better than the pentamer, and the dimers formed between two adjacent hexamers or between a pentamer and a hexamer can further lead to a more stable structure. These understandings could provide theoretical guidance for the modification of candidate vaccine with HBc-VLPs as the carrier.