In this dissertation, I combine physical modeling, mathematical analysis, and numerical computation to study several problems concerning electrostatic interactions in biological molecular systems. Biomolecules, such as DNAs and proteins, are highly charged. They interact with each other and with the mobile ions and the polarized solvent to generate strong forces that affect crucially many molecular processes, such as molecular conformational changes, recognition, and self-assembly. My goal is to develop rigorous theories and efficient computational tools to understand some of the principles and mechanisms underlying such complicated interactions. Specifically, I study three closely related problems. The first one is the effect of ionic sizes that has been experimentally observed to be significant. Using a size-modified mean- field model and extensive Monte Carlo simulations, I capture details of the competitive adsorption of counterions of multiple species to charged surfaces. The simulations confirm the crucial role of ionic valence-to- volume ratios that had been predicted by the mean-field theory. The second problem is how the ionic concentration dependent dielectric coefficient can affect the equilibrium properties of charged system. I construct a variational model, derive rigorously the first and second variations of the free-energy functional, and prove that the functional is nonconvex. Numerical studies reveal several important features that cannot be modeled with a constant dielectric coefficient. Finally, I study a general problem of the solvation of charged molecules, aiming at understanding how electrostatic interactions contribute to the shape of underlying molecules. I develop a phase-field implicit-solvent approach to minimize a free -energy functional that couples the interfacial energy, van der Waals solute-solvent interaction energy, and electrostatic energy. I also design numerical methods to implement this approach. A common theme of my dissertation 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