The present study applied atmospheric and room temperature plasma (ARTP) mutagenesis to improve Bacillus cereus strain for neutral proteinase production. The selected improved stable strain, Phe−Tyr-ARTP-60-E demonstrated a 2.88-fold enhanced proteinase yield and 2.81-fold enhanced proteinase activity with shortened fermentation time. An artificial neural network-embedded genetic algorithm (ANN-GA) was employed to optimize the enzyme's three-phase partitioning (TPP) recovery using Python software. The optimized solution of the hybrid model recommended 38.05% wv−1 (NH4)2SO4, 1.0:2.0 crude extract/tert-butanol ratio, pH 5.84 at 23.4 °C, for the exclusive partitioning of proteinase in the intermediate phase, with 202% recovery, 13-fold purity with specific activity of 922 Umg−1. Sensitivity analysis of the model revealed that enzyme recovery was most sensitive to crude extract/tert-butanol ratio with total sensitivity (ST) of 0.253 while purification factor was most sensitive to pH with an ST of 0.444. In silico analysis using Expasy ProtParam tool revealed that mutant proteinase had 3 amino acid substitutions; Val-to-Cys; Val-to-His and Ile-to-Lys, at positions 152, 288, and 292, respectively, and 39.2 kDa molecular weight, 6.42 isoelectric point, and 33.35 instability index. The 3-Dimensional structure of the enzyme predicted by ModWeb v r273 Server revealed 39% sequence identity with Pseudoalteromonas lipolytica thermolysin. The macromolecule increased transmittance in lysozyme-coated contact lenses by 99.6% in 45 min. We conclude that ARTP mutagenesis and ANN-GA optimized TPP are effective and efficient recalcitrant strain improvement and product recovery methods, respectively. Further genetic tools and analysis of the neutral proteinase-producing gene may be required for increased improvement toward sustainable application.