Abstract <p>An important side of biochemical reactions is the possibility of binding of several ligands to a biomolecular target. Multiple ligand attachment should be taken into account in the analysis of protein–ligand interactions and in evaluation of the drug distribution in living systems. This article focuses on the detailed molecular simulations of successive binding of two molecules of non-steroidal analgesic aspirin (As) at the known sites 1–3 of human serum albumin. The experimental data on multiple binding of aspirin to albumin are inconclusive. Docking of the As^– anion to albumin predicts that stability of the protein–ligand complexes changes in the order 1 &gt; 3 &gt; 2 for these sites. Subsequent molecular dynamics simulations have shown that the docked As^– positions at site 3 are not stable. The free energies of ligand binding Δ<i>G</i>_b have been calculated using linear interaction energy method with additional contributions of the ligand internal energy and the entropy of ligand binding. The calculations show that the most probable reaction path corresponds to binding of As^– at the Sudlow site 1 with Δ<i>G</i>_b1= ‒8.2 kcal&#xa0;mol^–1 and, after that, at the Sudlow site 2 with Δ<i>G</i>_b2= ‒4.5&#xa0;kcal&#xa0;mol^–1. The calculated values of Δ<i>G</i>_b agree with the known experimental data. The predicted stoichiometry of the albumin–As^– complexes is 2. Negative cooperative effect is found for binding of two As^– molecules with albumin. The employed computational approaches can be useful in molecular modeling of transport of multiple medicinal molecules by human serum albumin.</p>