Existing fundamental theorems for mean-square convergence of numerical methods for stochastic differential equations (SDEs) require globally or one-sided Lipschitz continuous coefficients, while strong convergence results under merely local Lipschitz conditions are largely restricted to Euler-Maruyama type methods. To address these limitations, we introduce a novel localized version of the fundamental mean-square convergence theorem for SDEs with locally Lipschitz coefficients, which naturally arise in a wide range of applications. Specifically, we show that if a numerical scheme is locally consistent in the mean-square sense of order $q>1$, then it is locally mean-square convergent with rate $q-1$. Building on this result, we further prove that global mean-square convergence follows, provided that both the exact solution and its numerical approximation admit bounded $2p$th moments for some $p>1$. These new convergence results are illustrated on a class of locally Lipschitz SDEs of Hodgkin-Huxley type, characterized by a conditionally linear drift structure. For these systems, we construct different Lie-Trotter and Strang splitting methods exploiting their conditional linearity. The proposed convergence framework is then applied to these schemes, requiring innovative proofs of local consistency and boundedness of moments. In addition, we establish key structure-preserving properties of the splitting methods, in particular state-space preservation and geometric ergodicity. Numerical experiments support the theoretical results and demonstrate that the proposed splitting schemes significantly outperform Euler-Maruyama type methods in preserving the qualitative features of the model.