We present the first application of deep learning for gravitational wave parameter estimation of binary black hole mergers evolving on quasi-circular orbits with aligned or anti-aligned spins. We use root-leaf structured networks to ensure that common physical features are shared across all parameters. In order to cover a broad range of astrophysically motivated scenarios, we use a training dataset with over $10^7$ modeled waveforms to ensure local time- and scale-invariance. The trained models are applied to estimate the astrophysical parameters of the existing catalog of detected binary black hole mergers, and their corresponding black hole remnants, including the final spin and the gravitational wave quasi-normal frequencies. Using a deterministic neural network model, we are able to efficiently provide point-parameter estimation results, along with statistical errors caused by the noise spectrum uncertainty. We also introduce the first application of Bayesian neural networks for gravitational wave parameter estimation of real astrophysical events. These probabilistic models were trained with over $10^7$ modeled waveforms and using 1024 nodes (65,536 core processors) on the Theta supercomputer at Argonne Leadership Computing Facility to reduce the training stage to just thirty minutes. In inference mode, both the deterministic and Bayesian neural networks estimate the astrophysical parameters of binary black hole mergers within 2 milliseconds using a single Tesla V100 GPU. Both deterministic and Bayesian neural networks produce agreeing parameter estimation results, which are also consistent with Bayesian analyses used to characterize the catalog of binary black hole mergers observed by the advanced LIGO and Virgo detectors.