Electrostatic interactions play an important role in many complex charged systems, such as biological molecules, soft matter material, nanofluids, and electrochemical devices. In this work, we develop mathematical theories and computational methods to understand such interactions, particularly in charged biological molecular systems. Our main contributions include: 1. Theoretical studies of mean-field variational models of ionic solution and that of molecular surfaces with the Poisson-Boltzmann electrostatics; 2. Design and implementation of the corresponding computational algorithms, and conduct extensive Monte Carlo simulations and numerical solutions of partial differential equations for charge-charge interactions; 3. Discovery of various interesting properties of charged molecules, validate some experimental results, and clarify some confusion in literature. A common theme of this work is the variational approach. Many physical effects such as ionic size effects, solvent entropy, concentration dependent dielectric response can be incorporated into a mean-field free-energy functional of ionic concentrations coupled with the Poisson equation for electrostatics. The techniques of analysis developed in this work may help improve the understanding of the underlying physical properties of charged systems and provide new ways of studying analytically and numerically other problems in the calculus of variations.