Towards A Novel Multi-Epitopes Chimeric Vaccine for Simulating Strong Immune Responses and Protection against Morganella morganii

Morganella morganii is one of the main etiological agents of hospital-acquired infections and no licensed vaccine is available against the pathogen. Herein, we designed a multi-epitope-based vaccine against M. morganii. Predicted proteins from fully sequenced genomes of the pathogen were subjected to a core sequences analysis, followed by the prioritization of non-redundant, host non-homologous and extracellular, outer membrane and periplasmic membrane virulent proteins as vaccine targets. Five proteins (TonB-dependent siderophore receptor, serralysin family metalloprotease, type 1 fimbrial protein, flagellar hook protein (FlgE), and pilus periplasmic chaperone) were shortlisted for the epitope prediction. The predicted epitopes were checked for antigenicity, toxicity, solubility, and binding affinity with the DRB*0101 allele. The selected epitopes were linked with each other through GPGPG linkers and were joined with the cholera toxin B subunit (CTBS) to boost immune responses. The tertiary structure of the vaccine was modeled and blindly docked with MHC-I, MHC-II, and Toll-like receptors 4 (TLR4). Molecular dynamic simulations of 250 nanoseconds affirmed that the designed vaccine showed stable conformation with the receptors. Further, intermolecular binding free energies demonstrated the domination of both the van der Waals and electrostatic energies. Overall, the results of the current study might help experimentalists to develop a novel vaccine against M. morganii.


Introduction
Antibiotic resistance (AR) is a global health crisis. AR is a subset of antimicrobial resistance (AMR) and happens when bacteria evolve mechanisms to withstand attacks by antibiotics. AR can evolve by natural courses forced by the continued misuse of antibiotics [1]. The resistance is mounting to seriously high levels across all countries of the world. Novel resistance mechanisms are evolving and spreading worldwide, making our efforts to treat common infectious diseases less effective [2]. As a consequence, infections caused by bacterial pathogens are becoming tough to treat, even sometimes impossible to 2 of 26 treat. According to the recent updates by the Center for Disease Control and Prevention (CDC), each year AR pathogens cause 2.8 million infections resulting in more than 35,000 deaths. In other words, it is one death every 15 min and one infection every 11 s [2][3][4]. In summary, we are losing the fight against bacterial pathogens. We are heading for a post-antibiotic era; therefore, urgent actions are needed to manage the AR crises. There are well established ways to lower the global burden of AR, such as the development of novel classes of antibiotics, the improvement of sanitation and hygiene, antibiotic stewardship, avoiding the routine use of antibiotics in agriculture and veterinary practice, and stopping its inappropriate use in treating viral infections [5].
Vaccination is an excellent alternative to combat AR bacterial pathogens. Historically, the use of vaccines as a tool to manage AR bacterial pathogens has been under-estimated, but its positive effect in reducing AMR is very well established [6]. As an example, Streptococcus pneumoniae (pneumococcal) conjugate and Haemophilus influenzae type B (Hib) vaccines have remarkable track records in reducing antibiotic use and AR as well as preventing life threatening diseases caused by these bacteria [7]. Therefore, new technologies for vaccine development will provide a potential solution to tackle AR. The growing amount of genomic data and advancements in bioinformatics tools have brought a revolution in the vaccine development process [8,9]. Reverse vaccinology (RV), a genome-based vaccine development pipeline, has contributed considerably to the identification of new vaccine candidates [8]. RV is pioneered by Dr. Rino Rappuoli and has revolutionized the vaccine development pipeline, particularly against pathogens for which Pasteur's principles of vaccinology have failed [10]. RV has been effectively used in meningococci serogroup B vaccine (4CMenB) development [11]. In the recent past, a recombinant chimeric peptide vaccine was designed in silico and evaluated experimentally to show 50% protection in hamsters against the infection [12]. Moreover, due to genomic diversity in bacterial pathogens, classical RV has been modified to pan-genomic based RV (PGRV) [13] to identify core proteome antigens. In the recent past, PGRV successfully mapped four protective antigens in Streptococcus agalactiae genomes [14].
This work involves core genomics, subtractive proteomics, and RV in combination with biophysical approaches to construct a multi-epitope vaccine and to decipher its binding potential with host immune system components as well as evaluate its potential in providing immune protection against M. morganii. M. morganii belongs to the family of Enterobacteriaceae and causes nosocomial infections, especially urinary tract and wound infections [15]. Some strains of the pathogen are resistant to oxacillin, penicillin, first-generation and second-generation cephalosporins, ampicillin/sulbactam, macrolides, fosfomycin, colistin, lincosamides, and polymyxin B [15]. In addition to this AR spectrum exhibited by M. morganii, there is no vaccine in development for the pathogen which may make the situation worse while treating these infections. Hence, substantial efforts are required to screen protective antigens from the pathogen genome that can be subjected easily to experimental evaluations. This in turn will save time and reduce the costs that usually go into the experimental vaccine candidate's prioritization. As AR is increasing in bacterial pathogens and there are hurdles in classical vaccinology, computer-aided vaccine design could provide easy access to surface-exposed protective antigens that otherwise use resources and are time consuming in experimental vaccine development. The study also employed B and T-cell epitope predictions, the analysis and processing of potential and safe antigens, population coverage and conservation analysis, toxicity prediction of the antigens, allergenicity evaluation, molecular docking, molecular dynamics simulation and binding energies estimations to evaluate the binding strength of the vaccine to host immune receptors. To accomplish the task of successful in silico vaccine design, several online web servers and bioinformatics tools were used.

Materials and Methods
The stepwise methodology followed for the design of a novel multi-antigenic epitope vaccine against M. morganii in this study is schematically shown in Figure 1.

Materials and Methods
The stepwise methodology followed for the design of a novel multi-antigenic epitope vaccine against M. morganii in this study is schematically shown in Figure 1.

M. Morganii Predicted Proteomes Retrieval
The predicted set of proteins of fully sequenced genomes (eight in number at time of the study) of M. morganii were retrieved from the NCBI [16] genome database and subjected to a bacterial pan-proteomic analysis [17][18][19] tool to extract the core proteins of the pathogen. Fast clustering of the proteomes was accomplished by setting the sequence identity cut-off at 50%, and the resulting file containing the core proteins was considered for further analysis as the proteins were conserved among all the strains.

Subtraction of Core Proteins
The subtractive proteomic approach was used for the analysis of the core proteins in order to identify potential vaccine candidates [20]. The subtractive proteomic approach is an in silico approach for the identification of vaccine targets which works by excluding all proteins which are not essential in the design of a vaccine candidate [21]. The first step in the subtractive proteomic approach was the removal of all paralogous proteins that were achieved through the Cd Hit analysis. To predict subcellular localization, non-homologues proteins were used in the PSORTb v3.0 [22] analysis. All virulent proteins were

M. Morganii Predicted Proteomes Retrieval
The predicted set of proteins of fully sequenced genomes (eight in number at time of the study) of M. morganii were retrieved from the NCBI [16] genome database and subjected to a bacterial pan-proteomic analysis [17][18][19] tool to extract the core proteins of the pathogen. Fast clustering of the proteomes was accomplished by setting the sequence identity cut-off at 50%, and the resulting file containing the core proteins was considered for further analysis as the proteins were conserved among all the strains.

Subtraction of Core Proteins
The subtractive proteomic approach was used for the analysis of the core proteins in order to identify potential vaccine candidates [20]. The subtractive proteomic approach is an in silico approach for the identification of vaccine targets which works by excluding all proteins which are not essential in the design of a vaccine candidate [21]. The first step in the subtractive proteomic approach was the removal of all paralogous proteins that were achieved through the Cd Hit analysis. To predict subcellular localization, non-homologues proteins were used in the PSORTb v3.0 [22] analysis. All virulent proteins were identified through BLASTp in the virulent factor database (VFDB) against a complete set of all the virulent proteins in the database [23]. The inclusion criteria were >30% sequence identity and a bit score of >100. Additionally, an antigenicity analysis was performed using the Vaxijen 2.0 online webserver [24]; the proteins defined as probably antigenic were those proteins whose antigenic score was greater than 0.5. The allergenicity of the proteins was predicted using Allertop 2.0 [25]. Next, all the non-redundant proteins were subjected to homology analysis against human and three Lactobacillus species: L. rhamnosus (taxid: 47715), L. casei (taxid: 1582), and L. johnsonii (taxid: 33959) with the selection criteria of a sequence identity of <30%, a bit score of > 100, and an E-value of 10 −4 to avoid auto immune reactions and the accidental inhibition of the good probiotic bacteria. This task was achieved via BLASTp [26]. Transmembrane helices were predicted using TMHMM 2.0 [27] with a cut off value of >1. Proteins with more than one transmembrane helix were discarded from further analysis [28]. Protein stability was examined using the Protparam tool [29].

Epitope Prediction Phase
In vaccine design, the prediction of B-cell and T-cell epitopes is crucial in order to elicit both cellular and humoral immune responses against the antigen [30,31]. B and T-cell epitopes were predicted for the shortlisted proteins from the previous step [32]. First, we predicted linear B-cell epitopes by using Bepipred linear epitopes 2.0 [33] and then B-cell epitopes were further used for the T-cell epitopes via the IEDB T-cell prediction tool [34] to determine the B-cell epitopes' binding ability for MHC class I and II. The method used for B-cell epitope prediction was IEDB recommended [35] while the peptides were selected on the basis of low percentile rank and if they were considered as a good binder to immune cell receptors [36]. Each selected epitope was analyzed for its binding affinity with DRB*0101 as this covers about 95% of the world population [37]. Antigenicity, allergenicity, toxicity, and solubility were checked for each of the selected epitopes by using the Vaxijen 2.0 [24], Allertop 2.0 [25], ToxinPred [38] and InvivoGen [39] tools, respectively. After all these analyses, the shortlisted epitopes were next subjected to multi-epitope vaccine design and processing.

Multi-Epitope Vaccine Design and Processing
One key issue in the single peptide base vaccine is that it cannot generate proper immune responses [40]. To overcome this problem, a multi-epitope-based vaccine was designed that comprised several different types of immunodominant epitopes rather than a single epitope [41]. A multi-epitope peptide vaccine construct containing different immunodominant epitopes is considered to be a good vaccine strategy to evoke substantial immune responses [42]. In order to make a multi-epitope construct, all the selected epitopes were linked with each other through a specific linker (GPGPG), and finally the designed multi-epitope vaccine construct was joined to the N-terminal of the beta subunit of the cholera toxin, which is considered a good and safe adjuvant [43].

Physiochemical Properties of the Final Vaccine Construct
The designed vaccine was further checked for physiochemical properties i.e., molecular weight, instability index, aliphatic index through an online Protparam (ExPasy) https://web.expasy.org/protparam/, (accessed on 10 May 2021) web server [29].

Structure Modelling of the Vaccine
The final multi-epitope vaccine construct was modeled ab initio for its 3D structure with the help of 3DPro [44]. Moreover, we re-validated the antigenicity and solubility of the vaccine using ANTIGENpro [44] and SOLpro solubility [44] using a Scratch protein predictor. Several loops of the vaccine were modeled through an online Galaxy WEB server http://galaxy.seoklab.org/, (accessed on 12 June 2021) [45]. After loop modelling, the loop modeled structure was submitted for refinement in the GalaxyRefine tool, which is available at http://galaxy.seoklab.org/, (accessed on 13 June 2021) [46]. The refinement was performed for structure errors and to lower the global binding energy of the vaccine.

Disulfide Engineering and Codon Optimization
The structural stability of the vaccine candidate was improved via disulfide engineering using Design 2.0 [47]. In disulfide engineering, the mutant structure was created by incorporating di-sulfide bonds in the vaccine structure. In order to obtain a high level of expression of the cloned vaccine sequence in the Escherichia coli system, the codon optimization approach was applied [48]. In this process, the sequence of the model vaccine construct was reverse translated to the DNA sequence through the Java Codon Adaptation Tool (JCat) [49].

Molecular Docking and Refinement
Molecular docking was performed using PATCHDOCK [50] and refined through the FIREDOCK server [51]. Molecular docking of the vaccine was performed with TLR4 (PDB: 4G8A), MHC-I (PDB ID: 1L1Y) and MHC-II (1KG0) receptors. In total, 20 docked solutions were predicted by PATCHDOCK that were ranked based on global binding energy. The FIREDOCK server re-ranked the solutions after removing many steric clashes and intermolecular conformational errors. The best conformation docked complex was visualized using UCSF Chimera 1.15 [52].

Molecular Dynamics Simulations
Molecular dynamic simulations of the docked complexes were performed for 250 ns to evaluate the structural stability of the systems. The simulations were carried out using the AMBER20 simulation package [53]. The simulation protocol described by the authors of [54] was followed to accomplish the assay. Briefly, the antechamber module [55] of AMBER was used to pre-process the systems while the parameters were defined using AMBER FF14SB force field [56]. The systems were solvated into TIP3P water box, where they were neutralized by adding appropriate amounts of counter ions. Afterward, the systems were subjected to energy minimization that can be split into hydrogen atoms energy minimization, water box energy minimization, and non-heavy atoms energy minimization. Next, the systems were gradually heated to 300 K and the temperature was maintained using Langevin dynamics [57]. Constraints on the hydrogen bonds were achieved using the SHAKE algorithm [58]. Equilibration of the systems was achieved for 100 ps, followed by pressure equilibration using NPT ensemble [59]. Last, the production run of simulations was performed for 250 ns on a time scale of 2 fs. The simulation trajectories were examined for different structural analysis using CPPTRAJ module [60] of AMBER. Visual inspection of the trajectories was completed using the Visual Molecular Dynamics (VMD) tool version 1.9.3 [61]. Each of the above steps was conducted with a different set of parameters as described by the authors in [54].

Binding Free Energies Calculation
The binding free energies of the docked complexes were calculated using the MM/PBSA and MM/GBSA approaches available in AMBER20 [62]. Both of these analyses were conducted using the MMPBSA.py module of AMBER [63]. Only 100 frames were considered while estimating the binding free energies.

C-immune Simulations
The immunogenic efficacy of the final vaccine construct was evaluated by performing in silico immune simulations with the help of C-immSim server 10.1 [64]. The server uses the position-specific score matrix (PSSM) and various other machine learning techniques to predict and study epitope and immune interactions.

Results
In this research work, a total of eight completely sequenced genomes of M. morganii were obtained from NCBI. Complete information about the strain's proteome can be found at the following link https://www.ncbi.nlm.nih.gov/genome/browse/#!/prokaryotes/10 874/, (accessed on 1 March 2021).

Bacterial Pan-genome Analysis (BPGA)
The bacteria strains had 16,880 core proteins collectively. The average number of core proteins encoded by each genome was around 2110. The number of accessory proteins, unique proteins, and exclusively absent proteins varied according to the strain, but the average values were 1164, 146, and 106, respectively. The number of proteins of each M. morganii strain is graphically presented in Figure 2A. The core-pan plot indicates that the predicted proteome of the pathogen is open and there is a high chance of gaining new genes over time due to genome plasticity. Moreover, COG distribution analysis reported that the core proteins were mostly engaged in metabolic biogenesis and regulation. The unique set of proteins (17,170 in number) were associated with the processing and storage of information. The information can be categorized into RNA processing, replication process, transcription and translation and recombination. Furthermore, the pan-phylogeny tree of the eight complete genomes of M. morganii was constructed, which is shown in Figure 2B. uses the position-specific score matrix (PSSM) and various other machine learning techniques to predict and study epitope and immune interactions.

Results
In this research work, a total of eight completely sequenced genomes of M. morganii were obtained from NCBI. Complete information about the strain's proteome can be found at the following link https://www.ncbi.nlm.nih.gov/genome/browse/#!/prokaryotes/10874/, (accessed on 1 March 2021).

Bacterial Pan-genome Analysis (BPGA)
The bacteria strains had 16,880 core proteins collectively. The average number of core proteins encoded by each genome was around 2110. The number of accessory proteins, unique proteins, and exclusively absent proteins varied according to the strain, but the average values were 1164, 146, and 106, respectively. The number of proteins of each M. morganii strain is graphically presented in Figure 2A. The core-pan plot indicates that the predicted proteome of the pathogen is open and there is a high chance of gaining new genes over time due to genome plasticity. Moreover, COG distribution analysis reported that the core proteins were mostly engaged in metabolic biogenesis and regulation. The unique set of proteins (17,170 in number) were associated with the processing and storage of information. The information can be categorized into RNA processing, replication process, transcription and translation and recombination. Furthermore, the pan-phylogeny tree of the eight complete genomes of M. morganii was constructed, which is shown in Figure 2B.

CD-Hit Analysis
The core proteins are a conserved set of the protein sequence which are shared by all the strains. The core proteins numbered 16,880 and were next analyzed for redundancy. This unveiled 14,924 redundant proteins and 1956 non-redundant proteins, as shown in Figure 3. The redundant sequences were duplicated and arose because of a duplication event during the evolution process. As such, for the computational vaccine design process these sequences were not required [65]. All the non-redundant proteins were subjected to subcellular localization and virulence analysis.

CD-Hit Analysis
The core proteins are a conserved set of the protein sequence which are shared by all the strains. The core proteins numbered 16,880 and were next analyzed for redundancy. This unveiled 14,924 redundant proteins and 1956 non-redundant proteins, as shown in Figure 3. The redundant sequences were duplicated and arose because of a duplication event during the evolution process. As such, for the computational vaccine design process these sequences were not required [65]. All the non-redundant proteins were subjected to subcellular localization and virulence analysis.

Subcellular Localization and Virulence Analysis
The human immune system can easily recognize those proteins which are present at the pathogen surface [66]. These surface proteins are also intrinsically more immunogenic and are in regular contact with the host cells. In 1956 core non-redundant protein sequences, 11 proteins were extracellular and outer membrane, and 14 proteins were found on the periplasmic membrane, as mentioned in Figure 3. Virulence refers to the relative degree of harmfulness of a pathogen to cause disease in a host [37]. After virulence analysis, among the core non-redundant proteins sequences, only 36 protein sequences were virulent in nature.

Antigenicity, Allergenicity, Human and Normal Microbiota Homology, and Transmembrane Helices and Stability Analysis
Antigenicity screening predicted 26 proteins as antigenic with a score of > 0.5 (Table  S1). To avoid allergic responses and auto immune reactions, all the virulent proteins were analyzed for allergenicity. The server found nine protein sequences as allergic among the total surface localized virulent proteins. Next, four proteins were human homologs and six proteins were identical to probiotic bacteria. The non-homology ensured that no unwanted auto-immune reactions will be generated if the proteins are used in vaccine design [28]. Similarly, the non-homology of the selected proteins also helps in avoiding the accidental inhibition of beneficial probiotic bacteria [67]. Transmembrane helices and the physiochemical analysis were checked and indicated that four proteins were removed from the study as they had transmembrane helices > 1. Further, eight proteins were predicted to be unstable (instability index is >40) with a molecular weight of >100 kDa as

Subcellular Localization and Virulence Analysis
The human immune system can easily recognize those proteins which are present at the pathogen surface [66]. These surface proteins are also intrinsically more immunogenic and are in regular contact with the host cells. In 1956 core non-redundant protein sequences, 11 proteins were extracellular and outer membrane, and 14 proteins were found on the periplasmic membrane, as mentioned in Figure 3. Virulence refers to the relative degree of harmfulness of a pathogen to cause disease in a host [37]. After virulence analysis, among the core non-redundant proteins sequences, only 36 protein sequences were virulent in nature.

Antigenicity, Allergenicity, Human and Normal Microbiota Homology, and Transmembrane Helices and Stability Analysis
Antigenicity screening predicted 26 proteins as antigenic with a score of > 0.5 (Table S1). To avoid allergic responses and auto immune reactions, all the virulent proteins were analyzed for allergenicity. The server found nine protein sequences as allergic among the total surface localized virulent proteins. Next, four proteins were human homologs and six proteins were identical to probiotic bacteria. The non-homology ensured that no unwanted auto-immune reactions will be generated if the proteins are used in vaccine design [28]. Similarly, the non-homology of the selected proteins also helps in avoiding the accidental inhibition of beneficial probiotic bacteria [67]. Transmembrane helices and the physiochemical analysis were checked and indicated that four proteins were removed from the study as they had transmembrane helices > 1. Further, eight proteins were predicted to be unstable (instability index is >40) with a molecular weight of >100 kDa as mentioned in the following Venn diagram in Figure 4. The shortlisted five proteins (TonB-dependent siderophore receptor) (serralysin family metalloprotease) (type 1 fimbrial protein) (flagellar hook protein (FlgE)) and (pilus periplasmic chaperone) which were non-allergic, nonhomologous to the host proteome, non-homologous to the probiotic bacteria, had no or <1 transmembrane helices and were within a range of molecular weights. These proteins are ideal candidates for subunit-based vaccine design. Further, the proteins were subjected to an epitope prediction phase in order to design a multi-epitope vaccine. mentioned in the following Venn diagram in Figure 4. The shortlisted five proteins (TonBdependent siderophore receptor) (serralysin family metalloprotease) (type 1 fimbrial protein) (flagellar hook protein (FlgE)) and (pilus periplasmic chaperone) which were nonallergic, non-homologous to the host proteome, non-homologous to the probiotic bacteria, had no or <1 transmembrane helices and were within a range of molecular weights. These proteins are ideal candidates for subunit-based vaccine design. Further, the proteins were subjected to an epitope prediction phase in order to design a multi-epitope vaccine.

B-Cell and T-Cell Epitopes Prediction
The active acquired immune responses are highly specific and specialized in clearing pathogens or inhibiting their growth [68]. Adaptive immunity basically generates memory B-cells that recognize the organism on successive encounters after initial recognition [69]. Such an immunological memory of adaptive immunity forms the foundation of vaccination. The B and T lymphocyte cells of the adaptive immunity are mainly involved in generating dependent antibodies and cellular immunity against invader organisms. Thus, in this study the final screened five protein sequences were used in B and Tcell epitope mapping.

B-Cell Epitope Prediction
The humoral immune response is referred to as the antibody-dependent immune response and it is activated when the B-cell matures and transforms into a plasma cell [70].

B-Cell and T-Cell Epitopes Prediction
The active acquired immune responses are highly specific and specialized in clearing pathogens or inhibiting their growth [68]. Adaptive immunity basically generates memory B-cells that recognize the organism on successive encounters after initial recognition [69]. Such an immunological memory of adaptive immunity forms the foundation of vaccination. The B and T lymphocyte cells of the adaptive immunity are mainly involved in generating dependent antibodies and cellular immunity against invader organisms. Thus, in this study the final screened five protein sequences were used in B and T-cell epitope mapping.

B-Cell Epitope Prediction
The humoral immune response is referred to as the antibody-dependent immune response and it is activated when the B-cell matures and transforms into a plasma cell [70]. The B-cell epitopes were predicted for one outer membrane (TonB-dependent siderophore receptor), three extracellular membrane proteins: serralysin family metalloprotease, type 1 fimbrial protein, flagellar hook protein (FlgE)), and one periplasmic membrane protein (pilus periplasmic chaperone). In total, 15 linear B-cell epitopes were predicted for TonBdependent siderophore receptor, 8 epitopes for serralysin family metalloprotease, 4 epitopes for type 1 fimbrial protein, 6 epitopes for flagellar hook protein (FlgE), and 2 epitopes were predicted for pilus periplasmic chaperone as tabulated in Table S2.

B-Cell-Derived T-Cell Epitope Prediction
Cellular immune responses, also known as T-cell-dependent immune responses, mainly function to kill infected cells [70]. T-cell lymphocytes in response to peptide antigens allow consequent multiplication and differentiation to constitute the primary immune responses [71]. To generate the cellular immune response, T-cell epitopes were predicted using B-cell epitopes called B-cell derived T-cell epitopes.  [71]. Each epitope was prioritized on the basis of low percentile rank. The lower the percentile score, the better the epitope is as a good binder, as tabulated in Table S3.

Epitope Prioritization Phase
In the epitope prioritization phase, all the predicted epitopes were next evaluated for DRB*0101 binding analysis, allergenicity, solubility and toxicity analysis. The binding affinity of the epitopes with the immune cell receptors is imperative; all the selected epitopes were checked for their ability to bind with the HLA DRB*0101 allele. This gene is prevalent in about 95% of the human population and any epitope binding to this allele will result in a better immune responses [37]. Only epitopes of IC50 values < 100 nM for DRB*0101 alleles were selected as they represent strong binding. The epitopes with IC50 values less than the threshold are shown in Table 1. Antigenicity was evaluated for each selected epitope with a cutoff value of 0.5, followed by allergenicity analysis. All probable antigenic and non-allergic epitopes were included in the study, while non-antigenic and allergic epitopes were excluded from study. Lastly, solubility and toxicity analyses were also performed in order to remove poor soluble and toxic epitopes to reduce solubility hurdles in the experimental evaluation of the vaccine as well as to avoid toxic resections.

Multi-Epitope Vaccine Construction
In a single peptide-based vaccine, one key issue is that it cannot generate proper immune responses. To overcome this problem, a multi-epitope-based vaccine was designed by linking different types of selected epitopes through specific GPGPG linkers [40]. Additionally, the epitope peptide was joined with the cholera toxin b subunit adjuvant with the help of an EAAAK linker. Both the used linkers are rigid and allow efficient separation of the epitopes so that they can be readily recognized by the host immune system. Similarly, the adjuvant used is safe and generates robust and specific immune responses against the antigen to which it is attached. The design of a multi-epitope vaccine is mentioned in Figure 5.

Structure Prediction, Loop Modelling, and Refinement
A 3D structure of the final vaccine construct was modeled using the sequence of the final model vaccine construct as is depicted in Figure 6. The structure modelling was completed ab initio rather than through homology or threading because no appropriate template structure was available. The vaccine was further deciphered as soluble and antigenic. Structural stability is one of the most important characteristics of a good vaccine candidate. Loop modelling was completed for the following residues: Met1-Thr11. Cys30-Ile38,Ser51-Gln70,Glu100-Lys112,Lys129-Thr148,Ser149-Gly169,Pro170-Gly185,Tyr189-Phe208,Gly209-Lys228,Gln229-Pro238,Gly242-Pro261,Gln262-Gly281,Pro282-

Structure Prediction, Loop Modelling, and Refinement
A 3D structure of the final vaccine construct was modeled using the sequence of the final model vaccine construct as is depicted in Figure 6. The structure modelling was completed ab initio rather than through homology or threading because no appropriate template structure was available. The vaccine was further deciphered as soluble and antigenic. Structural stability is one of the most important characteristics of a good vaccine candidate. Loop modelling was completed for the following residues: Met1-Thr11. Cys30-Ile38, Ser51-Gln70, Glu100-Lys112, Lys129-Thr148, Ser149-Gly169, Pro170-Gly185, Tyr189-Phe208, Gly209-Lys228, Gln229-Pro238, Gly242-Pro261, Gln262-Gly281, Pro282-Pro294, Gly304-Leu317, Gly321-Pro340, and Val341-Lys348. The refinement of the vaccine structure was completed to remove any high energy contact and to relax the structure. By doing so, we obtained five models as tabulated in Table S4. Model 1 was selected as it had better Rama-favored residue mapping (91.3%), an improved rotamer score (0.4), a clash score of 10.4, and a molProbity score of 2.04 compared to the initial structure. Pro294,Gly304-Leu317,Gly321-Pro340, and Val341-Lys348. The refinement of the vaccine structure was completed to remove any high energy contact and to relax the structure. By doing so, we obtained five models as tabulated in Table S4. Model 1 was selected as it had better Rama-favored residue mapping (91.3%), an improved rotamer score (0.4), a clash score of 10.4, and a molProbity score of 2.04 compared to the initial structure. Figure 6. Vaccine structure in 3D. The red color represents the adjuvant (cholera toxin B subunit), the purple color shows the vaccine construct, the green color is for the EAAAK linker, and the yellow color is for the GPGPG linkers.

Disulfide Engineering and Codon Optimization
In the design of a vaccine, improving the structural stability is an important objective. Disulfide bonds are covalent bonds that provide considerable structural stability to proteins and support the exact geometry of a given protein molecule [72]. In a multi-epitopebased vaccine, there is a chance that some of peptide residues will be enzyme degradable. To overcome this problem, disulfide engineering of the final vaccine construct was completed. In this step, enzyme degradable amino acid residues were replaced with cysteine residues, as in the mutant structure B shown in Figure 7. No significant changes were noticed in the structure of the vaccine after disulfide engineering. The VaxiJen antigenicity score of the disulfide engineered vaccine sequence was 0.8100. Moreover, all the epitopes were fully maintained in the disulfide engineering. Figure 6. Vaccine structure in 3D. The red color represents the adjuvant (cholera toxin B subunit), the purple color shows the vaccine construct, the green color is for the EAAAK linker, and the yellow color is for the GPGPG linkers.

Disulfide Engineering and Codon Optimization
In the design of a vaccine, improving the structural stability is an important objective. Disulfide bonds are covalent bonds that provide considerable structural stability to proteins and support the exact geometry of a given protein molecule [72]. In a multi-epitope-based vaccine, there is a chance that some of peptide residues will be enzyme degradable. To overcome this problem, disulfide engineering of the final vaccine construct was completed. In this step, enzyme degradable amino acid residues were replaced with cysteine residues, as in the mutant structure B shown in Figure 7. No significant changes were noticed in the structure of the vaccine after disulfide engineering. The VaxiJen antigenicity score of the disulfide engineered vaccine sequence was 0.8100. Moreover, all the epitopes were fully maintained in the disulfide engineering.
Codon optimization allows the improvement of gene expression and an enhancement of the translation efficiency by correcting the codon usage of a given sequence according to a host codon usage pattern [73]. The vaccine codon adaptation index (CAI) value was 0.92 and the GC content was 57.08%. Both these values are considered ideal for the expression process. Further, other factors were also evaluated and set to the non-binding site of the prokaryotic ribosome, the inactivation of rho-independent transcription termination, and the non-restriction enzymes cleavage sites. based vaccine, there is a chance that some of peptide residues will be enzyme degradab To overcome this problem, disulfide engineering of the final vaccine construct was co pleted. In this step, enzyme degradable amino acid residues were replaced with cyste residues, as in the mutant structure B shown in Figure 7. No significant changes w noticed in the structure of the vaccine after disulfide engineering. The VaxiJen antigenic score of the disulfide engineered vaccine sequence was 0.8100. Moreover, all the epito were fully maintained in the disulfide engineering.

Molecular Docking and Refinement
To generate a proper host immune response, the designed vaccine construct protein should interact with different types of the host's innate and adaptive immune cells and their receptors competently. A molecular docking analysis approach was applied to forecast the binding affinity of the designed vaccine construct with human immune cell receptors. A blind docking of the vaccine construct with TLR4, MHC-I, and MHC-II was performed [58]. First, the 3D structures of MHC-I, MHC-II, and TLR4 were retrieved and subjected to docking analysis. All of the interacting residues were assessed on the basis of the shape complementarity principle. This blind docking analysis is an essential step to evaluate the structure of a designed vaccine construct and to select all those epitopes with a high capability of interaction toward the selected immune receptors molecules. The results of the MHC-I, MHC-II, and TLR4 blind docking are shown in Table S5. In each case, 20 solutions were produced.
The PATCHDOCK results of the top 10 docked complexes were subjected to a refinement of the steric clashes. The complexes with the lowest global energy were ranked top, and selected further for binding mode and interaction studies through UCSF Chimera 1.13.1. For each receptor, the top docked solution was selected. In the case of MHC-I, solution nine was selected as it had the lowest global energy of −8.87 kJ·mol −1 with good contribution from attractive van der Waals (−20.43 kJ·mol −1 ), repulsive van der Waals (8.30 kJ·mol −1 ), atomic contact energy (ACE) (−3.02 kJ·mol −1 ) and hydrogen bond energy (−4 kJ·mol −1 ). Similarly, for MHC-II and TLR4, solution nine was selected based on the good global binding energy. The vaccine was observed to dock strongly to the receptors and formed strong intermolecular interactions. The FireDock refinement results for MHC-I, MHC-II, and TLR4 are provided in Table 2. The docked intermolecular conformation of the vaccine with MHC-I, MHC-II, and TLR4 is shown in Figure 8.

Chemical Interactions of the Vaccine with MHC-I, MHC-II, and TLR-4
Interaction between the vaccine and host immune cell receptors is very crucial to understand as it allows users to highlight the residues which are important in vaccine recognition and binding. The chemical interactions between the vaccine construct and TLR4, MHC-I, and MHC-II immune receptors were determined using the protein-peptide molecular docking approach, and the specific residues' interaction with MHC-I, MHC-II, and TLR 4 were checked in UCSF Chimera. The designed multi-epitope vaccine showed interactions with several residues of MHC-I within 3 Å. These interactions were both hydrophobic and hydrophilic. Similarly, the vaccine also produced a strong interaction network with the MHC-II molecule. All the interactions were within close distance and were of different types, including hydrogen bonding, salt-bridges, and van der Waals. TLR4 is one of the members of the TLR family and is a class of transmembrane peptides that basically belong to pattern recognition receptors (PRR) and are usually expressed on dendritic and macrophages cells. The interacting residues between the vaccine and MHC-I, MHC-II, and TLR4 are mentioned in Table 3.

Molecular Dynamics Simulation
Molecular dynamic simulation analysis basically checks the dynamic behavior of macromolecules [61]. The simulations analysis includes root mean square deviation (RMSD) [74], root mean square fluctuation (RMSF) [75], and radius of gyration (RoG) [76]. All these analyses were performed based on the carbon alpha atom of the complexes. These analyses were carried out to investigate whether the vaccine binding to the receptors was stable or not and whether the interactions remained intact throughout the simulation time. Stable binding of the vaccine with the receptors will ensure its proper presentation to the host immune cells that will further allow the activation of immune pathways to clear the antigen. The RMSD plot of the systems was uniform with no major structural changes observed. The TLR4-vaccine complex reported some deviations, but the systems achieved equilibrium towards the simulation end. The RMSD of the systems fluctuated around 5-6 Å ( Figure 9A). Second, the RMSF was calculated to shed light on the residue flexibility of the receptors in the presence of the vaccine molecule ( Figure 9B). The majority of the systems' residues were within a good stability (<3 Å). Some of the residues were found to show a higher flexibility, which is the outcome of the loop pressure on the system. However, these variations did not affect the vaccine binding to the receptors. Lastly, RoG analysis was performed to examine the systems' compactness verses time ( Figure 9C). As depicted in RMSD, all the systems were compact and did not experience any drastic variations.

Hydrogen Bonds Analysis
The H-Bonds plugin in VMD was used to count the number of hydrogen bonds formed throughout the simulation time as shown in Figure 10. All these bonds were within the cut-off distance of 3 Ǻ. A rich hydrogen bond clustering pattern was revealed between the vaccine construct and the TLR4, MHC-I, and MHC-II receptors. This depicts the strong formation of complex and high affinity of the vaccine to the receptors.

Hydrogen Bonds Analysis
The H-Bonds plugin in VMD was used to count the number of hydrogen bonds formed throughout the simulation time as shown in Figure 10. All these bonds were within the cut-off distance of 3Ǻ. A rich hydrogen bond clustering pattern was revealed between the vaccine construct and the TLR4, MHC-I, and MHC-II receptors. This depicts the strong formation of complex and high affinity of the vaccine to the receptors.

Determination of the Binding Free Energies
The binding free energies of the docked complexes were calculated using MM-GBSA and MM-PBSA approaches and were used to validate the binding ability of the docked complexes. The total binding free energies of the TLR4-vaccine complex, the MHC-I-vaccine complex, and the MHC-II-vaccine complex were −78.4 kcal/mol, −81.44 kcal/mol, and −55.01 kcal/mol, respectively. The major contributor to the net binding energy came from the van der Waals energy as well as electrostatic energy while the non-favorable contribution came from the solvation energy. The non-polar solvation binding energy also favored the complex formation between the vaccine and receptors. Details of the binding free energies for each complex are tabulated in Table 4.  ), EEL (electrostatic), EGB (polar solvation energy of MM-GBSA), ESURF (non-polar solvation energy), Delta G gas (net gas phase energy), Delta G solv (net solvation energy), Delta Total (net energy of system).

Determination of the Binding Free Energies
The binding free energies of the docked complexes were calculated using MM-GBSA and MM-PBSA approaches and were used to validate the binding ability of the docked complexes. The total binding free energies of the TLR4-vaccine complex, the MHC-Ivaccine complex, and the MHC-II -vaccine complex were −78.4 kcal/mol, −81.44 kcal/mol, and −55.01 kcal/mol, respectively. The major contributor to the net binding energy came from the van der Waals energy as well as electrostatic energy while the non-favorable contribution came from the solvation energy. The non-polar solvation binding energy also favored the complex formation between the vaccine and receptors. Details of the binding free energies for each complex are tabulated in Table 4. Table 4. MMGBSA/PBSA binding free energy results of the vaccine construct with MHC-I, MHC-II, and TLR4 complexes.  VDWAALS (van der Waals), EEL (electrostatic), EGB (polar solvation energy of MM-GBSA), ESURF (non-polar solvation energy), Delta G gas (net gas phase energy), Delta G solv (net solvation energy), Delta Total (net energy of system).

Computationally Immune Simulations
The immunogenic potency of the model vaccine construct was evaluated by performing a computational immune simulation with the help of the C-ImmSim server, which uses a position-specific score matrix (PSSM) and various other machine learning techniques to predict and study epitope and immune interactions. In the C-immune simulation analysis, upon the maximum level of vaccine antigen exposure to the human immune system for 50 days, an increase in the production of adaptive immune responses in the form of IgG and IgM antibodies was detected. The IgM antibody was also seen at a high level. Secondary immune responses followed by tertiary immune responses led to the maximum production level of B-cells and a high production of "IgM + IgG, IgM, IgG1 + IgG2, IgG1 and IgG2", as mentioned in Figure 11A. Likewise, the production of interferon gamma was greater than a 400,000 count per ml for almost 35 days as shown in Figure 11B. The different B and T-cell immune responses to the vaccine are presented in Figure S1. Similarly, the response of different immune cells to the vaccine is illustrated in Figure S2.

Computationally Immune Simulations
The immunogenic potency of the model vaccine construct was evaluated by performing a computational immune simulation with the help of the C-ImmSim server, which uses a position-specific score matrix (PSSM) and various other machine learning techniques to predict and study epitope and immune interactions. In the C-immune simulation analysis, upon the maximum level of vaccine antigen exposure to the human immune system for 50 days, an increase in the production of adaptive immune responses in the form of IgG and IgM antibodies was detected. The IgM antibody was also seen at a high level. Secondary immune responses followed by tertiary immune responses led to the maximum production level of B-cells and a high production of " IgM + IgG, IgM, IgG1 + IgG2, IgG1 and IgG2", as mentioned in Figure 11A. Likewise, the production of interferon gamma was greater than a 400,000 count per ml for almost 35 days as shown in Figure  11B. The different B and T-cell immune responses to the vaccine are presented in Figure  S1. Similarly, the response of different immune cells to the vaccine is illustrated in Figure  S2.

Discussion
The emergence of multi-drug resistant strains of bacterial pathogens offers a serious threat to human health [1,77]. This is particularly endangering the efficacy of the current list of antibiotics. This also impacts the development of new drugs, as the pharma industries have shown little interest because of challenging regulatory requirements and reduced economic incentives [78]. The AMR issue can be addressed by developing effective vaccines to stop bacterial infections [79]. Not all vaccinations are effective and helpful in eradicating targeted organisms, but they can at least reduce the infection burden that in turn will result in lowering antibiotic use and hence the evolution of novel resistant strains [79]. Traditional vaccinology approaches suffer from several shortcomings i.e., a lengthy processing time, the generation of inaccurate immune responses in the developed vaccine, high costs, reduced safety, less specificity, hypersensitivity, and less stability [80]. A substantial amount of genomic data is now available to help us in the identification of the antigenic proteins best suited for the development of novel vaccines [81]. The combination

Discussion
The emergence of multi-drug resistant strains of bacterial pathogens offers a serious threat to human health [1,77]. This is particularly endangering the efficacy of the current list of antibiotics. This also impacts the development of new drugs, as the pharma industries have shown little interest because of challenging regulatory requirements and reduced economic incentives [78]. The AMR issue can be addressed by developing effective vaccines to stop bacterial infections [79]. Not all vaccinations are effective and helpful in eradicating targeted organisms, but they can at least reduce the infection burden that in turn will result in lowering antibiotic use and hence the evolution of novel resistant strains [79]. Traditional vaccinology approaches suffer from several shortcomings i.e., a lengthy processing time, the generation of inaccurate immune responses in the developed vaccine, high costs, reduced safety, less specificity, hypersensitivity, and less stability [80]. A substantial amount of genomic data is now available to help us in the identification of the antigenic proteins best suited for the development of novel vaccines [81]. The combination of immunoinformatics and subtractive proteomics is a more attractive strategy in recent times to design low-cost and safe vaccines [82].
The current research study is based on an in silico approach for the design of a multiantigenic epitope vaccine against M. morganii. The study used an integrated approach including bacterial pan genome analysis, subtractive proteomics, epitope prediction and analysis, multi-epitope design and processing, receptor preparation, molecular docking, molecular dynamics simulation, C-immune simulations, and binding free energies calculations. In the initial steps, completely sequenced genomes of M. morganii were analyzed for core sequences and dispensable and unique sequences. Only core sequences were picked and processed as they represent proteins present in all the strains. This warranted the use of only high conserved sequences which would allow the development of a broad-spectrum vaccine [65]. In the subcellular localization analysis, outer membrane, extracellular membrane, and periplasmic membrane proteins were considered as vaccine targets [83]. Surface localized proteins are more potent for generating an immune response compared to other localized proteins. These proteins contain antigenic determinants and are in frequent contact with the host [84]. Moreover, they function as virulent proteins by activating immune signaling pathways [37]. Moreover, vaccine targets should be non-homologous to human and human intestinal microbiota to avoid autoimmune reactions and damage to the host beneficial bacteria [65]. In this study, only non-homologous, non-redundant, non-toxic, probable antigenic and immunogenic proteins sequences were subjected to epitope predictions. Antigenic and immunogenic proteins can provoke substantial and specific cellular and humoral immune responses. Five virulent and antigenic proteins (TonB-dependent siderophore receptor, serralysin family metalloprotease, type 1 fimbrial protein, flagellar hook protein, FlgE, and pilus periplasmic chaperone) were shortlisted as good vaccine candidates. Following that, an immunoinformatics pipeline was applied on the five vaccine targets to design a potential multi-epitope-based vaccine construct. After the epitope predictions, all the predicted epitopes were checked for further analysis including antigenicity, toxicity, adhesion probability, water solubility, and binding affinity with DRB*0101. In the multi-epitope vaccine construction phase, different types of epitopes were linked through specific GPGPG linkers because a single peptide-based vaccine cannot generate appropriate immune responses. Additionally, the designed epitope construct was joined to cholera toxin b subunit adjuvant with the help of an EAAAK linker in order to increase its antigenicity. A multi-epitope peptide plus an adjuvant allows robust B and T-cell responses [41]. The physiochemical properties of the vaccine were then checked to guide experimentalists in the experimental evaluation of the vaccine [85]. As the sequence homology template of the vaccine was not available, ab initio structure modeling of the vaccine was completed and refined for steric clashes so a proper conformation structure was used [86]. In the molecular docking phase, the designed vaccine construct was checked for its binding affinity with immune cell receptors [87] and validated by dynamics study [88]. A hydrogen bond analysis was also conducted to determine the intermolecular strength of the interactions between the vaccine and the receptors [89]. Hydrogen bonds are electrostatic forces formed when a hydrogen atom binds covalently to a more electronegative atom and is shared with another strongly electronegative atom. It is formed when hydrogen atoms are shared between electronegative hydrogen bond acceptors and donors. Mostly, hydrogen bonded to fluorine, oxygen, and nitrogen can either donate or accept hydrogen. For the binding free energies of the systems, the trajectories of the molecular dynamics simulations were subjected to MM/PBSA and MM/GBSA analysis [62]. All these predictions recommended that the designed vaccine could bind to the immune receptors efficiently and were capable of providing immune protection to the host. However, these findings must be evaluated experimentally in in vitro and in vivo models for biological potency against M. morganii. A previous study conducted by the authors of [54] on the design of an in silico multi-epitope vaccine evaluated the vaccine potency using different computer-aided vaccine strategies and revealed that the vaccine could generate strong immune responses against invader Enterobacteriaceae. Another similar approach designed by the authors of [90] on the development of a chimeric vaccine against Proteus mirabilis found that broad-spectrum multi-antigenic, non-redundant, conserved, and surface localized peptides can provoke specific and accurate immune responses. While designing a vaccine, they analyzed three different types of proteins, (AtfC, PMI2533, and PMI1466). All these proteins fulfilled promising vaccine parameters and were able to serve good subunit vaccine candidates. Moreover, another study used a computer-aided vaccine design strategy against Pseudomonas aeruginosa [91], Providencia rettgeri [92], Streptococcus pneumoniae [11], and Klebsiella pneumoniae [93]. Computational vaccine designing strategies are rapidly emerging mainly because of the exponential growth of genomic data. Such analyses are highly specific, effective, and can guide the production of a safe vaccine against a number of microbial pathogens [93][94][95]. The vaccine design herein fulfilled all the parameters of a good vaccine. However, the study had limitations in relation to experimental testing which needs to be performed in the future to confirm vaccine immune protection potency against M. morganii.

Conclusions and Limitations
In this study, we used different applications of computer-aided vaccine design including RV, subtractive proteomics, immunoinformatics, and different biophysical analysis to propose a novel multi-epitope vaccine to train the human immune system to fight against M. morganii, which is a nosocomial bacterial pathogen responsible for several infections. The pathogen is evolving new resistant mechanisms, and no licensed vaccine is available to curtail its infections. Our vaccine construct was comprised of multi-epitopes picked from potential vaccine proteins that were prioritized based on several studies in the literature that have reported vaccine properties. The designed vaccine consisted of antigenic and non-toxic virulent epitopes predicted using core, non-redundant, host non-homologous, antigenic, and experimentally favored proteins. The fusing of the epitopes and the adjuvating was completed using specific GPGPG and EAAAK linkers. The designed vaccine showed excellent binding to the immune receptors, revealing fittest binding conformation, and generating robust binding energies. We believe that the findings and predictions of the current study will speed up the vaccine development process against M. morganii and will deliver data that might fast track vaccine development against said organism. Furthermore, the outcomes of the study will also be able to save time and millions of dollars, and the in silico designed vaccine construct would be helpful for experimental vaccinologists in the formulation of a vaccine against M. morganii infections in both prophylactic and therapeutically circumstances. Though we remained very strict regarding the selection criteria at each step of the study, some limitations still need to be overcome in future studies. First, the order of the epitopes in the vaccine construct is something that needs thorough experimental evaluation to obtain the best combination for the maximum level of immune response. Second, the MHC epitope prediction algorithms are under the process of refinement and as such we were not sure about those predictions. Lastly, the real immune protection of the vaccine required extensive in vivo and in vitro testing. Moreover, the study has some recommendations that need to be considered in future computational vaccine design strategies. Rapid developments have been observed in recent times across the biological and computational sciences that need to be properly integrated to further advance computational vaccine design. Although advancements have been reported in genome sequencing and development of novel bioinformatics tools for vaccine design, more efforts are still needed to improve the epitope prediction. The implementation of artificial intelligence and machine learning in computational vaccine design is needed. As bacterial diversity is very high across the globe, the availability of all prevalent strains of bacterial species must be included in the vaccine design process to ensure the design of a broad-spectrum vaccine rather than strain-specific vaccine. Only very limited sequence data on AR bacterial pathogens from Pakistan are available in international databases. This must be highlighted to obtain complete sequence bacterial genomes from our country that can be utilized in vaccine design.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/ijerph182010961/s1. Figure S1.  Table S1. Antigenic exo-proteome, secretome, and periplasmic proteins. Only proteins with a cut-off value of ≥0.5 were selected as antigenic and presented in the table. Table S2. B-cell epitopes predicted from antigenic proteins. The epitopes were variable in length and have bepipred linear epitope prediction scores of ≥0.50. Epitopes with a residue length of <9 were not selected for further analysis. Table S3. MHC-I and MHC-II epitopes predicted from B-cell epitopes. A percentile rank is the result of comparing the given peptide's predicted binding affinity against a set of similarly sized peptides selected randomly from the SWISS-PROT database. All the presented epitopes had the lowest percentile score for a given B-cell epitope (cut-off < 100). Epitopes that are common in MHC-I and MHC-II were selected for further processing. Table S4. Refinement models of the loop model vaccine structure. The model selection was based on the overall lowest GDT-HA (global distance test-high accuracy) score, lowest RMSD (root mean square deviation) score, lowest MolProbity score, lowest clash score, lowerst number of poor rotamers, and highest number of residues in the Ramachandran-favored regions of the Ramachandran plot. For each parameter, there was not specific threshold value. Table S5. Docking score of the top 20 complexes of the designed vaccine construct to MHC-I/MHC-II/TLR-4 generated by PATCHDOCK server. The complexes were ranked on PATCHDOCK score. A higher score implies the best docked complex and vice versa. Energy is presented in kJ·mol −1 (lower energy solution is often regarded as the best docked complex). ACE stands for atomic energy contact.