<p>The nonlinear Poisson–Boltzmann equation (NPBE) is an elliptic partial differential equation used in applications such as protein interactions and biophysical chemistry (among many others). It describes the nonlinear electrostatic potential of charged bodies submerged in an ionic solution. The kinetic presence of the solvent molecules introduces randomness to the shape of a protein, and thus a more accurate model that incorporates these random perturbations of the domain is analyzed to compute the statistics of quantities of interest of the solution. When the parameterization of the random perturbations is high-dimensional, this calculation is intractable as it is subject to the curse of dimensionality. However, if the solution of the NPBE varies analytically with respect to the random parameters, the problem becomes amenable to techniques such as sparse grids and deep neural networks. In this paper, we show analyticity of the solution of the NPBE with respect to analytic perturbations of the domain by using the analytic implicit function theorem and the domain mapping method. Previous works have shown analyticity of solutions to linear elliptic equations with interfaces but not for nonlinear problems. We further show how to derive a priori bounds on the size of the region of analyticity. This method is applied to the Cucurbita Maxima Trypsin Inhibitor I (CMTI-I) molecule to demonstrate that the convergence rates of the quantity of interest are consistent with the analyticity result. Furthermore, the approach developed here is general enough to be applied to other nonlinear problems in uncertainty quantification.</p>