A methodology for large scale molecular dynamics simulation of a solvated polarizable protein, using a combination of permanent and inducible point dipoles with fluctuating and fixed charges, is discussed and applied to the simulation of water solvated bovine pancreatic trypsin inhibitor (BPTI). The electrostatic forces are evaluated using a generalized form of the P3M Ewald method which includes point dipoles in addition to point charge sites. The electrostatic configuration is propagated along with the nuclei during the course of the simulation using an extended Lagrangian formalism. For the system size studied, 20000 atoms, this method gives only a marginal computational overhead relative to nonpolarizable potential models (1.23-1.45) per time step of simulation. The models employ a newly developed polarizable dipole force field for the protein(1) with two commonly used water models TIP4P-FQ and RPOL. Performed at constant energy and constant volume (NVE) using the velocity Verlet algorithm, the simulations show excellent energy conservation and run stably for their 2 ns duration. To characterize the accuracy of the solvation models the protein structure is analyzed. The simulated structures remain within 1 Å of the experimental crystal structure for the duration of the simulation in line with the nonpolarizable OPLS-AA model.