Identification of potential mutational hotspots in serratiopeptidase to address its poor pH tolerance issue

Abstract Serratiopeptidase is the multifunctionality metalloendopeptidase extensively employed in biopharmaceutical and industrial biotechnology. Despite its poor pH tolerance, serratiopeptidase must withstand the highly acidic environment of the gastrointestinal tract to be used as a potent anti-inflammatory and analgesic medication. In earlier studies, post-translational deamination related mutations showed alteration in the net charge of protein’s surface. Therefore, the current study aimed to enhance the acid resistance of serratiopeptidase via implementing computational interventions to screen out the most stable mutational hotspot. The methodology used in this study is as follows: (a) Higher accessibility to surface (b) 4 Å away from active site region to avoid interference with its proteolytic activity, and (c) By converting non-conserved amide residues to acidic residues. A docking study has been conducted to establish the substrate specificity and binding affinity to native and mutant proteins. The docking outcomes were then validated using molecular dynamic simulations to clarify each mutant’s molecular stability and conformation while preserving their activity. The results showed that N412D is the best-screened mutant with negative electrostatic potential that can alter the overall charge on the protein’s surface with increased H+ ions. Alteration in overall charge leads to protein surface more acidic that causes a common ion effect in stomach pH and act as a buffer which could stabilize the serratiopeptidase amid extreme pH. Communicated by Ramaswamy H. Sarma


Introduction
Biotechnology has broadened the applicability of traditional enzymes, allowing them to be employed in industries as diverse as food, leather processing, detergent, pharmaceutical, and many others (Fasim et al., 2021).Recent advancements in protein/enzyme engineering and site-directed progression have empowered us to modify enzymes with the novel trend or even for unique processing parameters (Zhu et al., 2011).Proteases account for 20% of the 60% of enzymes sold worldwide.Serratiopeptidase is a microbial metalloprotease that is employed therapeutically and industrially.This peptidase was first discovered in the 1960s from Enterobacterium Serratia E-15, an isolate from the intestine of the silkworm Bombyx mori.In silkworms, serratiopeptidase aids in the disintegration of the cocoon, allowing the silk moth to emerge (Miyata et al., 1970).Serratiopeptidase is also referred as serrapeptase, serratia-protease, and serralysin.It is a metalloendoprotease that cleaves hydrolytic bonds far from the substrate's termini.
Serratia marcescens, a facultative rod-shaped gram-negative anaerobe, emits serrapeptase, an alkaline protease.This protease exhibits two domains architecture and a zinc (Zn 2þ ) metal ion in its catalytic site, which is localized in the N-terminal (proteolytic domain) region attached to the three histidines from the HEXXHXXGXXH motif (His176, His180, and His186), along with glutamate (Glu177) a catalytic base and tyrosine (Tyr216).The protein's C-terminal domain is made up of repeat-in toxin, a glycine-aspartate-rich sequence that is responsible for the binding of seven calcium ions to the protein (Hamada et al., 1996;Molla et al., 1986) (Figure 1).It has been demonstrated that it shows higher activity at pH 9.0 with a temperature of 40 � C and gets deactivated within 15 minutes at 55 � C (Santhosh Kumar, 2018).The gene encoding serratiopeptidase is composed of 471 amino acids (Serra471), having a molecular mass ranging from 40 to 60 kDa (Indrayati et al., 2014;Kaviyarasi et al., 2016;Nakahama et al., 1986).Serratiopeptidase is an extensively used serine protease with several beneficial applications, including anti-inflammatory, anti-oedema, anti-bacterial, analgesic, immunomodulatory, antiviral, and fibrinolytic properties.It has been scientifically tested to digest arterial plaques, cysts, and blood clots, as well as dissolve dead tissues without disrupting the live tissue surrounding the inflamed area (Ethiraj & Gopinath, 2017;Isono et al., 1974;Mazzone et al., 1990;Yamasaki et al., 1967).Serratiopeptidase lowers vascular permeability mediated by histamine, bradykinin, and serotonin; degrades pathological secretions of proteins, and promotes the absorption of destroyed products via plasma and lymphatic fluid (Nakamuha et al., 2003).With the growing demand for serratiopeptidase in the medicinal and industrial domains, we must work to reduce all restraints such as ADR (adverse drug reactions), pH instability in acidic environment of the stomach, and low permeability through the intestinal mucosa.In previous literature, Selan and colleagues performed mutagenesis to evaluate the anti-biofilm property of serratiopeptidase in the absence of glutamic acid by substituting it with alanine (Selan et al., 2015).In an another study to enhance the thermal stability and protease activity of serratiopeptidase four truncated mutants and one disulfide bridge was generated (Rouhani et al., 2021).So far as we know, no scientific study has addressed the low pH tolerance issue of such an extensive natural enzyme.Serratiopeptidase is stable in alkaline pH but degrades in the acidic pH of the stomach, limiting its use as a potent therapeutic agent.(Rouhani et al., 2021).
In this study, we addressed the pH tolerance issue of serratiopeptidase by using in-silico mutagenesis to screen stable mutations to enhance serratiopeptidase's role as a potent therapeutic medication without eliciting any side effects like NSAIDs.In general, in-silico methodologies appear to be an efficient way to carry out a full evaluation of the anticipated consequence of mutations (da Silva J� unior et al., 2021;Pinto et al., 2019).To gain a better understanding of how mutations affect the stability of wild protein, we used a variety of computational tools, with a focus on covering the key factors driving protein stability.The goal of this work is to verify computationally that among all mutants, which mutations have proven to be stable across various servers of mutation prediction.Our findings based on molecular docking and dynamics will provide the additional details for exploring these mutations in wet-lab studies.

Dataset preparation
As of now, four crystal structures of Serratia protease have been identified.The 1SAT (Baumann, 1994), 1SMP (Baumann et al., 1995), 1SRP (Hamada et al., 1996), and 5D7W (Wu et al., 2016) from Serratia species have been submitted to the protein database (PDB) with resolutions of 1.8, 2.3, 2.0, and 1.1 Å, respectively.Among four PDB structures, 1SAT is the only structure with no mutation that possesses 1.8 Å resolution.We have chosen 1SAT as our template framework for the further steps.We mutated the (amide residues) asparagine (Asn) and glutamine (Gln) residues on the surface of 1SAT to (acidic residues) aspartic (Asp) and glutamic (Glu) residues.A protein's net charge can be altered by substituting amide residues with acidic ones.Thus, it affects the acid resistance issue on the protein's surface.Using the PyMOL (http://www.pymol.org)program, the identified amino acid substitutions were individually integrated into the existing 3 D enzyme model.The residues were chosen based on parameters such as non-conserved and not close to the active site (within a 4 Å radius).

In-silico screening of potential mutational sites
To see the stability of generated mutants compared to wildtype (WT) structure, we used a variety of published in-silico mutation prediction algorithms and web servers to forecast stability changes caused by single-point substitutions, primarily utilizing machine-learning methodologies.These web servers could be used to find out point mutations for stabilizing impacts to modify enzymatic stability rationally (bioinformatics, 2020;Cheng et al., 2006;Rodrigues et al., 2018).In this study, 6 predictive computer programs were used to predict stable mutations: Hot-spot wizard (Sumbalova et al., 2018), DynaMut server (Rodrigues et al., 2018), CUPSAT (Parthiban et al., 2006), I-Mutant 2.0 (Capriotti et al., 2005), ProtParam server (Gasteiger et al., 2005), and SWISS-MODEL (Waterhouse et al., 2018).Each computational method performed using default parameters to classify variants, as provided in the supplementary section (Tables S1-S6).

Surface accessibility percentage
In this study, we targeted the amide residues Asn and Gln on the protein surface for mutagenesis.Furthermore, posttranslational deamination of Asn and Gln residues can potentially change proteins net charge (Bischoff & Schl€ uter, 2012).To improve the acid resistance of serrapeptase, we altered the Asn and Gln residues on the surface of Serra471 to Asp and Glu, respectively (Liu et al., 2019).The HotSpot Wizard tool has been used to determine the surface accessibility of these residues.This tool generates detailed primary structure assessments and supports protein researchers in the rational design of site-directed mutations and targeted library resources.After entering the query structure, hotspot uses nineteen prediction algorithms and three databases for protein annotation (Sumbalova et al., 2018).In this process, acidic residues (Asp þ Glu) with surface accessibility of 50% or higher have been selected and altered with amide residues (Asn þ Gln) with the same attributes.These generated mutants would alter the charge on the protein's surface, which changes the pKa value by causing protonation and deprotonation, which can modify the pH value of the system.If the pKa value is lesser than the pH value, it undergoes deprotonation or vice-versa (Reijenga et al., 2013).

Missense mutation prediction
We used the DynaMut server to evaluate the effects of missense mutations on protein mobility, stability, and dynamics.DynaMut combines two unique well-established methodologies in which Normal Mode Analysis (NMA) illustrates protein conformational dynamics.The other approach, to see the effect of mutation on proteins stability and function (Bauer et al., 2019).In 1SAT, the overall stability prediction outcome was anticipated by applying the server's consensus approach for mutation prediction by aggregating findings from three servers, namely mutation Cutoff Scanning Matrix (mCSM), site-directed mutagenesis (SDM), and DUET server.This server computes the observed inconstancy in the protein structure (kcal mol À 1 ) and the anticipated variations in thermodynamic entropy between the WT and mutant forms (kcal mol À 1 K À 1 ).The DDG � 0 shows a stabilizing impact, whereas variants with DDG < 0 promote protein instability (Rodrigues et al., 2018).

Amino-acid substitution stability
We utilized one more program Cologne University Protein Stability Analysis Tool (CUPSAT) to assess protein stability upon substitutions.This program also indicates DDG, the variance in free energy generation between WT and mutant proteins.Output holds details about the structural properties (secondary structure, solvent accessibility, and torsion angles), its substitution site, and specific information on protein stability differences for the possible 19 amino acid mutations (Parthiban et al., 2006).

Prediction of protein stability
To predict protein stability changes upon different mutations, I-Mutant 2.0 (Capriotti et al., 2005), a support vector machine (SVM) based server, was utilized.This tool estimates modifications in protein stability induced by single-point substitutions.In the current study, we employed the protein sequence of 1SAT to make the prediction.This server was used as a predictor to denote the magnitude of a protein stability shift caused by mutation and as a logistic indicator to ascertain the DDG (relative change in Gibbs-free energy) (Wang et al., 2020).

Physicochemical characterization
The ExPASy ProtParam tool has been employed to examine the serratiopeptidase's physical and chemical properties (Gasteiger et al., 2005).This prediction server has been used to compute the physicochemical characteristics of the WT and mutated variants, such as molecular weight, amino-acid composition, theoretical isoelectric point (pI), atomic constitution, extinction coefficient, projected half-life, instability index, number of charged residues, aliphatic index (AI), and grand average of hydropathy (GRAVY).

3D Structure generation
The protein structure has been modelled, both WT and mutant, and a Ramachandran plot was being constructed using the homology-modelling database SWISS-MODEL (Waterhouse et al., 2018).This server enables the comparative analysis of the query/mutant and the protein's existing 3 D structures.

Protein-peptide docking
After estimating mutational hotspots, we docked wild and mutant-type enzymes with their substrate.Docked structures of WT and its mutants with substrate bradykinin (peptide) were constructed using the High Ambiguity Driven Docking (HADDOCK) program (De Vries et al., 2010).This server uses biochemical or biophysical interaction data to drive molecular docking.During this program, we predict active site residues of both enzyme and its substrate using Consensus Prediction Of interface Residues in Transient complexes (CPORT) web-server (de Vries & Bonvin, 2011).We have uploaded proteins (wild and mutants) and bradykinin in PDB format in the input dialogue box.PDB structures with ions such as zinc and calcium as their cofactors have given a charge and a HETATM (hetero-atom) flag.The docking method in HADDOCK is governed by Ambiguous Interaction Restraints (AIRs), which is essentially a set of algorithms that operate in collaboration with Crystallography and NMR System (CNS) and Ambiguous Restraints for Iterative Assignment (ARIA).AIRs derive from the most recent experimental data on the residues participating in the interaction between the two subunits.
In the distance restraints method, a fraction of the ambiguous restraints were randomly excluded with 2.0000 no partitions.The HADDOCK docking process begins with energy minimization and rigid body orientation randomization.We computed 1000 complex systems for simulated annealing and selected the 200 solutions with the least intermolecular interactions.After that, the optimized structures were improved in explicit solvent.For its scoring function, HADDOCK uses a three-step protocol: (a) Rigid body docking and randomization of orientation of complex featuring restraint-guided energy reduction method (it0) Quasi refinement in torsion angle space with backbone atoms and flexible side-chains (it1) Enhancement of the explicit water (itw) At the end of each phase, findings were graded using the weighted sum of thermodynamic properties that constitute the HADDOCK Score.This docked score is a summation of terms whose weights vary according to the different phases of the docking protocol: where E vdw and E elec depict non- , and E AIR is the restraint violation energy displaying the alliance between empirical and back-calculated statistical information.To further generate schematic diagrams of protein-peptide interactions, Ligplot tool was used (Wallace et al., 1995).

Molecular dynamic simulations
All-atom simulations were carried out for native and designed structures using the GROMACS v.4.6.7 package on a Linux kernel with GROMOS96 43a1 as the force-field parameter (van Gunsteren et al., 1996).The MD procedure has been initiated using the crystal structure of serratiopeptidase and generated mutants.Each protein system has placed in an SPC water box, and the radius from the water box's surface to all of the solute's atoms was set to 1.0 nm (Berendsen et al., 1981).
The system was subsequently neutralized by adding an appropriate number of counter ions (Na þ and Cl -) to balance the net charges of the system.The steepest descent algorithm was used to minimize the energy of the systems.After energy minimization, equilibrated the solvent and ions around the protein.
Equilibration is often conducted in two phases.In Phase-I system was equilibrated under an isothermal-isochoric or canonical (NVT) ensemble for 2 ns time steps.In output control, attribute data from coordinates and log files are saved and updated every 0.2 ps.In the next phase, the equilibration step was conducted under the isothermal-isobaric (NPT) ensemble.Again, this system was subjected to a 2 ns MD simulation using the leap-frog integrator with an input time step of 0.002 ps.The temperature of the system was set to 300 K by using velocity rescaling algorithm (Berendsen thermostat) with a thermostat time constant of 0.1 ps (Berendsen et al., 1984).The coupling method (Parrinello-Rahman pressure) with a 2 ps barostat time constant was employed to keep the system's pressure at a constant 1 bar (Parrinello & Rahman, 1981).The simulations were carried out for 100 ns, and the atomic coordinates were saved every 2 ps.MD trajectories analysis and visualization were conducted using the XMGRACE and visual molecular dynamics program (Humphrey et al., 1996).The resulting trajectory analysis was performed in terms of root mean square deviation (RMSD), root mean square fluctuation (RMSF), solvent accessible surface area (SASA), hydrogen bonds (H-bonds), and the secondary structure of both WT and mutants computed using gmx command do_dssp, where dictionary of secondary structure of proteins (DSSP) is the database of secondary structure predictions.

Electrostatic charge potential
Advanced Poisson-Boltzmann solver (APBS) v1.5 has been used to determine the hypothetical electrostatic potential of computation in wild and mutant (N412D) (Baker et al., 2001).
We utilized the PDB2PQR method to add missing side-chain atoms, hydrogens and assign radii and charges to individual atoms to prepare the molecule (Dolinsky et al., 2007).The grid spacing for map calculation was set to 0.50, with a surface visibility range of þ/À 5.00, omitting the solvent surface.In the command line parameter, the AMBERff99 force field was used (Spasic et al., 2012).The electrostatic potential visualization has been carried out using PyMOL v2.5.2 with different colour annotations red for negative charge, white for neutral, and blue for positive charge ranges between À 5 and þ5.

Screening of potential hotspots in the dataset using insilico prediction algorithms
In this research, comprehensive in-silico analysis has been conducted utilizing several computational tools comprising software and servers.The structure of WT and generated mutants were thoroughly investigated the mutations to draw a conclusion.Serrapeptase comprising of 54 acidic residues (Asp (39) þ Glu (15) ¼54), and 58 amide residues (Asn (35) þ Gln (23) ¼ 58).To accomplish post-translation deamination, a HotSpot Wizard tool was used to choose variants with point mutations on the surface of protein.These surface residues were important for the protein to show the area of protein which will directly in contact with the environment of gastrointestinal tract.
We have selected optimal mutational hotspots, which include nine Asp residues (N20, N25, N191, N262, N264, N289, N300, N346, N412) and four Gln residues (Q23, Q144, Q146, Q396) based on their secondary structure, accessible surface area (%), mutability score, relative flexibility, and B-factor.In the supplementary section, we have constructed a tabular dataset to illustrate all of these parameters (Table S1a and b).Missense mutation predictions have been made using data from the DynaMut server.This server generates results using a consensus strategy that incorporates three servers (mCSM, SDM, and DUET) and normal mode analysis (NMA).In the outcome of this server, three Q to E mutations (Q144E, Q146E, and Q396E) and seven N to D mutations (N20D, N191D, N262D, N264D, N300D, N346D, and N412D) were chosen in total (Table S2a and b).Additionally, we used CUPSAT server to determine the overall stability shifts and probable DDG (kcal/mol) values for all feasible substitutions.The server's major predictions revealed that the mutants Q144E and Q396E demonstrated stable mutation with an unfavorable torsion angle; conversely, predictions for Q146E were not stable and favorable.In contrast, N264D, N300D, and N346D demonstrated significant stability outcomes with a favorable torsion angle.Additionally, it was projected that the torsion percentages of N20D, N262D, and N412D were favorable, whereas N191D was predicted to be a stable but unfavorable torsion percentage.As a result, only two mutants with Q to E mutations and seven with N to D mutations have been chosen for further predictions (Table S3a and b).Furthermore, I-Mutant2.0was used to forecast the value of the (DG) free energy stability shift caused by a single point mutation in protein structure.The parameters for the prediction of each mutant protein stability variations were determined by setting the temperature at 40 � C and pH at 1.5 (extremely acidic) (Miyata et al., 1970, Table S4a and b).In the subset of more reliable outcomes sorted out based on the prediction mentioned above, only one mutant (Q396E) exhibited a substantial drop in stability, whereas the other eight mutants (one(Q-E) and seven (N-D)) exhibited a significant increase.

Physicochemical characterization and 3-D model prediction
The primary structure characterization of native and mutant structures has been estimated using the ExPASy tool ProtParam and tabulated in (Table S5a and b).Here, physicochemical features such as total no. of atoms, molecular weight (MW), isoelectric point (pI), negatively and positively charged residues, extinction coefficient, instability index, aliphatic index, and Grand Average hydropathy (GRAVY) indices of all the eight mutants were analyzed and compared with the native one.The MW of WT was 50275.50Da, theoretical isoelectric point 4.54, and an aliphatic index of 68.22.The instability index (II) was 22.54, and the hydropathy index was À 0.367.As we compared these parameters with the mutant structures, we found almost comparable MW that was 50275.50 for wild and 50276.48 for all the eight screened mutant forms.The theoretical pI value for native protein was less than 7, indicating the more significant number of positively charged amino acids.While in the case of mutants, we have observed an increase in negatively charged residues due to mutated asp and glu (-vely charged) residues.Furthermore, we calculated the extinction coefficients for all mutant proteins at 280 nm wavelength, demonstrating that all mutants were discovered to be stable.We have also evaluated the instability index of the selected mutants in comparison to WT.If the instability index is above 40 the mutants were considered to be unstable.Whereas, in our case all the mutants exhibited stable instability index ranges between 22.36 and 23.21.Further, all mutants aliphatic index (68.22),including the native one, demonstrated that they were stable across a wide temperature range.After physicochemical characterization, all the eight mutant structures were modelled using SWISS-MODEL server.The WT protein (1SAT) used as a template for all of the mutations, with 99.79% sequence identity.Further, we used Ramachandran plot to get insight into structural modifications of the predicted protein structures fitting in the favored or unfavored region.The native protein had a MolProbity score of 0.94, 97.65% of the residues were in the allowed area, 0% in the Ramachandran outliers section, with a clash score of 0.44% (Table S6a and b).All the selected mutants were in a favorable zone, comparable MolProbity score, with a stable Z-score not greater than À 4.0, indicating a poor quality model.To improve the dependability of quality estimates, the global Model Quality Estimate (GMQE) server used which provides overall model quality measurements ranging from 0 to 1.The value of model quality measurement has been reported to be less than one (�0.97) in eight mutants.In addition, we have also validated the modelled structure with artificial intelligence driven server AlphaFold v2.1.0.The result of superimposition showed 0.51 Å RMSD between the compared structures modelled by SWISS-MODEL and AlphaFold servers.The low RMSD values reaffirmed that both modelled structures did not render any disparity (Figure S1).Table 1 summarizes the overall predictions of all servers used, along with their outcomes.

Protein-peptide binding mode and interaction analysis
As stated in the introduction, serratiopeptidase also serves as a fibrinolytic enzyme, capable of digesting insoluble proteins such as fibrin without damaging healthy tissues.Furthermore, as per the literature, Serrapeptidase catalyzes the hydrolysis of inflammatory mediators (bradykinin, histamine, and serotonin), which aids in the pain reduction and oedema as well as the improvement of the vascular system, which in turn supports and significantly contributes to proteolytic, analgesic, and fibrinolytic activities (Desser et al., 1993).Serrapeptase acts exclusively on the Gly-Phe linkage in bradykinin peptide, exhibiting substrate specificity (Miyata et al., 1970).We used HADDOCK server for protein-peptide molecular docking.A molecular docking study is an optimization method that identifies the binding energy, intermolecular interactions, and inhibitory potential of a biological target of interest (Costa et al., 2020;Rego et al., 2022).Initially, we docked bradykinin to nine protein structures including WT and eight mutant forms.Prior knowledge of the active site of the serratiopeptidase allows us to predict where and how bradykinin will bind to the native surface.The H176, E177, H180, H186, and Y216 residues were selected as the active site residues attached to Zn 2þ (cofactor in native) using CPORT server (Figure 1) (Srivastava et al., 2017).The overall HADDOCK score has significant role in assessing the interaction between protein and peptide.The higher the value of the HADDOCK score, the greater the affinity between protein-peptide complexes (Van Dijk et al., 2006).Statistically, HADDOCK generated several clusters, but we have selected the structure with the best HADDOCK score for further analysis.
The HADDOCK score was nearly same in all the mutant complexes except for N191D which showing lower binding energy.The reason for lower binding energy was the 3.6 Å shift in the position of N191D to form an H-bond with the ligand's P3 residue, which was responsible for closing the opening to the active site cleft (Baumann, 1994) In addition, the RMSD value for all mutant complexes represents the mean deviations from the native one (Table 2).Other major contributing energy components for complexes are van der Waal energy and binding energy.Our data represented that van der Waal energy were higher (-ve) for all the complexes except Q144E and N300D, possessing À 37.2 þ/À 1.0 and À 33.7 þ/À 4.2 kcal/mol, respectively, where more negative van der Waal energy showed suitable complex formations.To test the binding energy between complexes, we choose the optimal structure that adopts the ideal pose for our purpose.We have selected only one complex with the lowest binding energy score from the top cluster comprising of four structure complexes.All mutants possessed a substantial binding energy score except for Q144E and N20D (À 157.31 and À 196.523 kcal/mol).The other parameter, buried surface area (BSA), represents the protein's surface that has not been exposed to solvent.A higher BSA score means a more compact and rigid biomolecular structure.However, in our scenario, our methodology requires more surface accessibility of residues to modify the total charge and pH of the original protein.Across all the mutant complexes, N262D and N412D (1100.0þ/À46.7 and 1169.4 þ/À 19.8) showed the least buried surface area compared to the wild one, i.e. 1175.1þ/À26.4.An additional attribute Z-score was examined, a statistical assessment where a cluster sits with the average score of all clusters.The set with the top position out of four will always have the lowest Z-score.All the complexes exhibited appropriate Z-score values as compared to the WT.
Furthermore, to study the protein-peptide interaction, Ligplot was used to generate the schematic 2-D representations, which facilitates the rapid inspection of intramolecular interactions and their strengths, including hydrophobic interactions and H-bonds.In contrast to WT complex (6 H-bonds), N20D and Q144E mutant complexes showed only 5 and 2 Hbonds respectively Figure S2.Apart from N20D and Q144E, all the other mutants were selected for further investigations, because their H-bonds and binding energy were higher or equal to the native complex.According to docking data, all mutant complexes exhibited favorable outcomes, as shown in Table 2, except for N20D and Q144E, which had the most negligible interactions.Further, based on molecular docking attributes, we have selected the six best mutant complexes (N191D, N262D, N264D, N300D, N346D, N412D) with WT to examine their conformational stability using MD analysis Figure 2.

Structural stability analysis
We performed 100 ns MD simulations on the six selected mutants and WT protein obtained from docking analysis to study their conformational space and dynamic processes of structural alterations.To estimate mutant's conformational fluctuations and backbone flexibility, we used open source software GROMACS to perform simulations (Medeiros et al., 2019;Neto et al., 2022).The RMSD measures the spatial variation among protein's initial and simulated structures.RMSD value for the c-A backbone was calculated to check the stability of the selected structures (Bhardwaj & Purohit, 2021;Kumar et al., 2022).All the selected mutants and WT structures showed average RMSD values between 0.30 to 0.37 nm (Figures S3 and 3A).Following that, all the mutants and WT proteins displayed equivalent divergence, with the backbone RMSD of �0.25 to �0.37 nm in the time frame of 15 ns to 75 ns, except for N300D and N346D mutants.From 75 ns onwards, the mutant showed slight fluctuations by achieving a backbone RMSD of �0.3 to �0.39 nm and then attained a stable trajectory similar to the WT.The RMSD analysis inferred that the generated trajectories were stable with a slight deviation in the starting of the MD simulations.Thus, providing an appropriate framework for additional evaluation.
Additionally, we assessed the RMSF of backbone c-A atoms to examine perturbations and backbone mobility in protein at the residue level (Figure 3B).The lower RMSF value implies less drift, whereas the larger RMSF value indicates flexibility (Kumar Bhardwaj et al., 2022;Singh et al., 2022b).The trajectories of wild and all the mutant types were studied, and they demonstrated that the key residues 176-186 (green dotted box) of Region-I showed lower fluctuations which is a significant factor for protein-peptide interactions.The RMSF values for all the trajectories were in the range of 0.05-0.1 nm.Furthermore, the Region-II Van der Waals energy À 39.9 þ/À 3.1 (peptidase domain) containing mutated residues showed RMSF values between 0.08-0.24nm.All the mutants exhibited almost identical trajectories to the WT with RMSF values ranging from �0.05 to �0.61 nm.All the mutated residues (acidic in nature) were existed in the peptidase domain and showed higher accessibility (%) to the protein's surface (Table S1(a) and Figure S4).Higher accessibility of the acidic residues on the surface of the protein is responsible for the reduction of positive charge.It protects protein from deteriorating at higher acidic pH (highly-vely charged residues) of the stomach by exhibiting a common ion effect.(Liu et al., 2019).However, aiming to enhance pH tolerance of serratiopeptidase, the mutant N412D (purple dotted box) shows higher accessibility (88%) compared to other mutants (Table S1(a)) with stable trajectories to the backbone of C-A atoms.The outcomes of RMSF analysis suggested that the acquired flexibility in peptidase domain residues as a result of N412D mutation was the cause of a significant docking score with an overall stable trajectory (Table 2 and Figure 3B).In addition, to check the robustness of the selected protocol we have performed the MD simulations in triplicate (Figure S5).The average deviation ranges between 0.2 to 0.3 nm in all the three simulations run and that is why we have selected the first trajectory.

Surface accessibility surface area and intermolecular H-bond analysis
SASA of the protein was determined and tracked throughout simulations to ensure the changes in the accessibility of protein within its solvent molecule.We calculated SASA and performed a comparative analysis to see the structural deviations of native and mutant structures.SASA is a crucial component in solubility and folding of protein.It is differentiated by a putative centre of the liquid sphere that interacts with the molecule's interface via van der Waals interactions (Gromiha et al., 2018).The surface area around the protein was determined for 100 ns, and it was observed that all the mutants had a SASA value that was fairly comparable to the WT, with subtle variations.All the mutant forms except N412D showed lesser accessibility to the solvent molecule  than the native structure, as shown in Figures S6 and 3C.
The SASA value of N412D complements the accessible surface area value obtained from the HotSpot Wizard server, as shown in Table S1.Further, we elaborated the finding of N412D mutant in contrast to WT.The initial structure of both proteins near 5 ns remains unchanged, with SASA area ranging from 185 to 197 nm 2 .In the native form, it drops steeply in the range 178 to 187 nm 2 across a time frame of 8 to 30 ns.At the same time, the mutant exhibited an overall increase in the SASA area and moved parallel to the x-axis (Figure 3C).Although the previous analysis revealed that the mutants had comparable trajectories to the WT, we also examined the mutants intermolecular H-bonds in comparison to the native one.H-bond analysis is critical in molecular recognition because it validates the functionality, stability, and specificity of interactions.Except for N412D, all mutants had average H-bonds in the range of 361-363 compared to WT, with �371 (H-bonds) (Figures S7 and 3D).The H-bonds in the N412D were comparable to those in the native.The analysis of intermolecular H-bonds revealed that the N412D mutant had a more stable conformation and fewer perturbations than the other mutant forms.

Secondary structure prediction using DSSP
We computed and compared the secondary structure of both N412D and WT proteins by executing the DSSP algorithm for each time frame.For secondary structure identification, DSSP employs the geometric (C-coil, b-sheets, bend, turn, and A-helices) and H-bond pattern based on electrostatic interactions (Singh et al., 2022a).The stability of the secondary structure was determined in the graph over a time interval of 100 ns.Interestingly, there was a rise in b-sheets in the mutant type from 35 to 45 ns and from residue 37 to 50 through 10 to 100 ns time frame.Our mutation region was showed in the graph between 400 and 450.
There was an increase in b-sheets in the mutant protein after changing Asn to Asp residue at position 412 (Figure 4).As the number of b-sheets increased, so did the number of Hbonding patterns in the peptide backbone, which further strengthened the secondary structure.The findings demonstrated a surge in H-bond turn situated primarily on the surface of the protein and possessing -ve charged residues with stable secondary structure, which satisfied the criteria of our proposed study.The outcomes of the DSSP investigation demonstrated that N412D had fewer perturbations and a stable trajectory with adequate secondary structure than WT.

Electrostatic charge potential analysis
Using the PyMOL programme, we noticed that the charge distribution becomes more negative with each time step, inferring that the number of H þ ions on the surface of the mutant protein rises in contrast to the WT.The presence of more H þ ions (protonation) signifies that the protein's surface is becoming more acidic, and the global region of the protein surface is also influenced, with charges shifting (Ward et al., 2013).In comparison to WT, all the snapshots of N412D mutant demonstrated desirable outcomes at different time steps ranging from 0 to 100 ns.On the other hand, the snapshots at the time intervals of 10 ns, 30 ns, 60 ns, and 70 ns produced excellent findings by displaying more negative regions near the mutational site.The -ve electrostatic potential near the mutational hotspot increased rapidly at 10 ns (red).Globally, it is nearly identical to the native protein, with the exception of the residues 106-161 becoming more negative.At 30 ns, the charge density on the WT protein shifted from neutral to more negative, this alteration changes the net charge of the surface.Furthermore, at around 60 ns, the protein's surface along the mutant region turned neutral (white) to negative (red).At 70 ns, the protein surface nearest the mutation revealed a strongly negative surface by reducing the neutral charge.Generally, it affects the region from residue 90 to 116 by shifting the charge from positive (blue) to negative (red) (Figure 5).Electrostatic charge analysis revealed that the N412D mutation influences the overall structure of a protein by increasing the -ve charge potential to the surface, which could act as a protective layer for serratiopeptidase degradation in the stomach's highly acidic environment (Figure S8).Such more mutations on the surface of protein could ameliorate the degradation issues of any therapeutically essential protein in the acidic conditions of gastrointestinal tract.

Conclusion
The present manuscript attempts to validate computationally stable mutational hotspots in serratiopeptidase.The current study's findings pave the way for rapid screening of robust serratiopeptidase variants for acid resistance in the gastrointestinal mucosa.According to the mutational tools and MD analysis outcomes, none of the six designed mutants differed significantly from the native protein.However, when all the mutants trajectory was analyzed and compared, N412D revealed the highest resemblance to the WT with adequate pattern of secondary structure, increased (-ve) electrostatic potential, and overall stable conformations.This research will unfold the pH-related issues that would enable us to economically and industrially endorse the usage of microbial peptidases as biopharmaceuticals around the world.

Figure 1 .
Figure1.A structural representation of the active site residues (His176, Glu177, His180, His186, and Tyr216) involved in interaction with the catalytic metal ion (zinc) and water molecule as equatorial ligand (aqua blue).

Figure 2 .
Figure 2. Graphical representation of protein-peptide interactions: A.) Figure depicting a complex of native serratiopeptidase protein (cyan) and its substrate bradykinin in stick (red) radiating a plot indicating the intramolecular interactions, and B.) Two dimensional interactions between substrate bradykinin and serratiopeptidase mutants.

Figure 3 .
Figure 3. Molecular dynamics simulations of WT (black) and N412D (red): A.) RMSD plots for backbone c-a atoms relative to their initial minimized structure, B.) RMSF plot for backbone c-a atoms relative to their initial minimized structure as a function of the residues, C.) Total solvent accessible surface area of the protein, and D.) Number of H-bonds in N412D and in WT protein.

Figure 4 .
Figure 4. Assessment of secondary structure elements of N412D in contrast to WT protein using dictionary of secondary structure of proteins analysis.

Figure 5 .
Figure 5. Electrostatic potential analysis of WT and N412D mutant at different time step using advanced Poisson-Boltzmann solver of PyMOL.

Table 1 .
Tabular representation of data using computer based algorithm for prediction of protein stability.

Table 2 .
Statistical study of protein-peptide docking analysis using HADDOCK server.