In this dissertation, we mainly discuss the image approximation methods to the reaction field and their applications to electrostatic interactions in molecular dynamics simulations. The Poisson-Boltzman (PB) equation, a three-dimensional second order nonlinear elliptic partial differential equation arising in biophysics, as well as the Debye-Huckel theory are considered as fundamentals throughout this work.We first outline a fourth-order image approximation proposed by Deng and Cai to the reaction field for a charge inside a dielectric sphere immersed in a solvent of low ionic strength [36]. To present such a reaction field, the image approximations employ a point charge at the classical Kelvin image charge location and two line charges that extend from the Kelvin image charge along the radial direction to infinity. A sixth-order image approximation is then developed, using the same point charge with three different line charges. Procedures on how to discretize the line charges by point image charges and how to implement the resulting point image approximation in O(N) complexity for potential and force field calculations are included. Numerical results demonstrate the sixth-order convergence rate of the image approximation and the O(N) complexity of the fast implementation of the point image approximation.We then apply the image-based reaction field method to the calculation for electrostatic interactions in molecular dynamics simulations. To extend a model developed by Lin et al. [30], a new hybrid solvation model, termed the Image-Charge Solvation Model (ICSM), is extended for simulations of biomolecules in ionic solvent, which combines the strengths of explicit and implicit solvent representations. In our model, an accurate and efficient multiple-image charge method to compute reaction fields is employed together with the fast multipole method for the direct Coulomb interactions. To minimize the surface effect, we use the periodic boundary condition strategy for nonelectrostatic interactions. We test our model in a simulation of sodium-chloride-water solvent. Using the Particle Mesh Ewald (PME) simulations as a reference, our results demonstrate that the proposed model can faithfully reproduce known solvation properties of sodium and chloride ions as well as many structural and dynamic properties of the water. We conclude that the employed model achieves convergence and controlled accuracy with only one image charge in the case of ionic solvent.