Molecular dynamics simulations have been used to investigate the thermodynamic stability of axial contacts in sickle-cell hemoglobin (HbS). Free energy changes were evaluated for the point mutation beta 121 Glu --> Gln in the axial contact region of HbS crystals. The calculations predict a free energy change of-3.6 kcal/mol per contact for the mutation, which is in qualitative agreement with experimental observations of aggravated sickling found in the double mutant Hb D Los Angeles (beta 6 Glu --> Val. beta 121 Glu --> Gln) relative to HbS (beta 6 Glu --> Val). The beta 121 Glu is sequestered in a salt link with beta 17 Lys located on the same polypeptide chain, making the Glu interactions with its surroundings similar in aggregates and individual hemoglobins. Due to this cancellation of the large electrostatic Glu contributions, the weak nonspecific interactions between the Gln and the neighboring polypeptide chain are the main contributing factor to the enhanced aggregation of Hb D Los Angeles relative to HbS. Together with the previous study of the lateral contact [K. Kuczera et al. (1990) Proceedings of the National Academy of Science USA, Vol. 87, pp, 8481-8485], the present results provide a more complete picture of the forces driving the sickling aggregation. A comparison of different treatments of internal flexibility in free energy simulations and analysis of rate of convergence of the different calculated properties has also been performed.