We have carried out a physicochemical and computational analysis on the stability of the intercalated structures formed by cytosine-rich DNA strands. In the computational study, the electrostatic energy components have been calculated using a Poisson-Boltzmann model, and the non-polar energy components have been computed with a van der Waals function and/or a term dependent on the solvent-accessible surface area of the molecules. The results have been compared with those obtained for Watson-Crick duplexes and with thermodynamic data derived from UV experiments. We have found that intercalated DNA is mainly stabilized by very favorable electrostatic interactions between hydrogen-bonded protonated and neutral cytosines, and by non-polar forces including the hydrophobic effect and enhanced van der Waals contacts. Cytosine protonation electrostatically promotes the association of DNA strands into a tetrameric structure. The electrostatic interactions between stacked C.C+ pairs are strongly attenuated by the reaction field of the solvent, and are modulated by a complex interplay of geometric and protonation factors. The forces stabilizing intercalated DNA must offset an entropic penalty due to the uptake of protons for cytosine protonation, at neutral pH, and also the electrostatic contribution to the solvation free energy. The latter energy component is less favorable for protonated DNA due to the partial neutralization of the negative charge of the molecule, and probably affects other protonated DNA and RNA structures such as C+-containing triplexes.