Glycosaminoglycans (GAGs) are anionic, periodic, linear polysaccharides which are composed of periodic disaccharide units. They play a vital role in many biological processes ongoing in the extracellular matrix. In terms of computational approaches, GAGs are very challenging molecules due to their high flexibility, periodicity, predominantly electrostatic-driven nature of interactions with their protein counterparts and potential multipose binding. Furthermore, the molecular mechanisms underlying GAG-mediated interactions are not fully known yet, and experimental techniques alone are not always sufficient to gain insights into them. The aim of this study was to characterize protein-ion-GAG complexes for the systems where ions are directly involved in GAG binding. Molecular docking, molecular dynamics and free energy calculation approaches were applied to model and rigorously analyse the interactions between annexins (II and V), calcium ions (Ca2+) and heparin (HP). The computational data were examined and discussed in the context of the structural data previously reported by the crystallographic studies. The computational results confirm that the presence of Ca2+ has a tremendous impact on the annexin-HP binding site. This study provides a general computational pipeline to discover the complexity of protein-GAG interactions and helps to understand the role of ions involved at the atomic level. The limitations of the applied protocols are described and discussed pointing at the challenges persisting in the state-of-the-art in silico tools to study protein-ion-GAG systems.