Cell membranes are composed of lipid bilayers that behave like a fluid in the in-plane direction---allowing molecules to rearrange without any static resistance---and behave like a solid in the out-of-plane direction---offering static bending resistance. This allows them to take a myriad of morphological shapes depending on their function. There are several crucial cellular processes like clathrin-mediated endocytosis, vesicle budding, and organelle fission, which rely on cell membranes undergoing large curvature driven deformations like budding, necking and tethering. As the cell membranes are usually only a few nanometres thick but can span a few hundred micrometers in their overall size, they can be effectively modeled as continuum surfaces. A finite element formulation with the Helfrich bending model as the strain energy functional is a common setup for such studies. However, the classical Helfrich model exhibits instabilities when modeling large deformations involving saddle-like geometries (negative Gaussian curvature), such as the catenoid shapes characteristic of membrane necks during fusion and fission. Additionally, as a consequence of the Gauss-Bonnet theorem, the Gaussian curvature is considered to be an invariant and neglected in simulations without a change in the surface topology or boundary. To address these limitations, we propose a modified Helfrich model that introduces a robust dependence on the Gaussian curvature with linear scaling just like in the classical model. We present a rigorous comparative analysis of the classical Helfrich energy, existing deviatoric curvature models, and our proposed formulation through both analytical and numerical benchmarks. The numerical framework is established using Isogeometric Analysis (IGA) based on Kirchhoff-Love shell theory, which naturally satisfies the C1-continuity requirements of the fourth-order bending problem while accurately capturing complex, non-axisymmetric deformations. Results demonstrate that the modified model successfully stabilizes intermediate saddle geometries and accurately captures energy barriers in large-deformation pathways where the classical model fails.