A general method has been developed to include the electrostatic reaction field in Brownian dynamics (BD) simulations of ions diffusing through complex molecular channels of arbitrary geometry. Assuming that the solvent is represented as a featureless continuum dielectric medium, a multipolar basis-set expansion is developed to express the reaction field Green's function. A reaction field matrix, which provides the coupling between generalized multipoles, is calculated only once and stored before the BD simulations. The electrostatic energy and forces are calculated at each time step by updating the generalized multipole moments. The method is closely related to the generalized solvent boundary potential [Im , J. Chem. Phys. 114, 2924 (2001)] which was recently developed to include the influence of distant atoms on a small region part of a large macromolecular system in molecular dynamics simulations. It is shown that the basis-set expansion is accurate and computationally inexpensive for three simple models such as a spherical ionic system, an impermeable membrane system, and a cylindrical pore system as well as a realistic system such as OmpF porin with all atomic details. The influence of the static field and the reaction field on the ion distribution and conductance in the OmpF channel is studied and discussed.