Molecular Dynamics Simulations in Designing DARPins as Phosphorylation-Specific Protein Binders of ERK2

Extracellular signal-regulated kinases 1 and 2 (ERK1/2) play key roles in promoting cell survival and proliferation through the phosphorylation of various substrates. Remarkable antitumour activity is found in many inhibitors that act upstream of the ERK pathway. However, drug-resistant tumour cells invariably emerge after their use due to the reactivation of ERK1/2 signalling. ERK1/2 inhibitors have shown clinical efficacy as a therapeutic strategy for the treatment of tumours with mitogen-activated protein kinase (MAPK) upstream target mutations. These inhibitors may be used as a possible strategy to overcome acquired resistance to MAPK inhibitors. Here, we report a class of repeat proteins—designed ankyrin repeat protein (DARPin) macromolecules targeting ERK2 as inhibitors. The structural basis of ERK2–DARPin interactions based on molecular dynamics (MD) simulations was studied. The information was then used to predict stabilizing mutations employing a web-based algorithm, MAESTRO. To evaluate whether these design strategies were successfully deployed, we performed all-atom, explicit-solvent molecular dynamics (MD) simulations. Two mutations, Ala → Asp and Ser → Leu, were found to perform better than the original sequence (DARPin E40) based on the associated energy and key residues involved in protein-protein interaction. MD simulations and analysis of the data obtained on these mutations supported our predictions.


Introduction
Protein kinases play a principal regulatory role in nearly all aspects of cell biology. The human genome encodes 538 protein kinases [1]. Out of all the post-translational modifications (PTMs), protein phosphorylation is the most widespread class used in signal transduction. One of the members of the protein kinase family, mitogen-activated protein kinase (MAPK) is a type of protein kinase that is specific to the amino acids serine and threonine (i.e., a serine/threonine-specific protein kinase). MAPK is involved in the most fundamental pathway to cell biology, known as the MAPK pathway, and it plays a crucial role in integrating cell surface signals to transcriptional regulation of the proteome [2]. The MAPK pathway is also referred to as the RAS-RAF-MEK-ERK signal cascade. The main function of this cascade is to regulate physiological processes by transmitting upstream signals to its downstream effectors. They are mainly involved in cell differentiation, proliferation, survival, and death. Targeting the MAPK pathway is thought to be a hopeful strategy for cancer therapy as it is the most often mutated signalling pathway in human cancer. In the past decades, extensive efforts have been made by different research groups, leading to the clinical success of BRAF and MEK inhibitors. In the pharmaceutical industry, computational drug design has played a vital role in the discovery, design, and analysis of drugs. Computer technology nowadays is so rich and advanced that the accuracy of biomolecular simulations is consistently high enough to be used to truly drive preclinical drug discovery projects. Among the computational tools used for drug discovery, molecular dynamics simulations (MDS) and related methods are routinely used nowadays. Their main contribution is towards the understanding of structural flexibility together with entropic effects of complex systems. This allows an in-depth estimation of the thermodynamics and kinetics related to drug-target recognition and binding [18]. MD simulations have been extensively used in the study of proteinprotein and protein-ligand interactions, and the study of the mechanism of drug action [19][20][21][22][23]. The plethora of applications of molecular dynamics simulations extend from the study of complex and dynamic processes that play a central role in biological systems to the structure determination from X-ray crystallography and NMR experiments. In biological systems, the main application of MDS focuses on studies related to the stability of proteins through conformational changes and folding. MD studies also enable the molecular recognition of cellular components such as DNA, membranes along with complexes (drug-receptor), and ion transport [24][25][26][27]. MD studies have also enabled the study of the mechanism of drug resistance [28][29][30][31]. Molecular dynamics simulations calculate the binding free energy, which is very helpful in investigating receptor-ligand interactions [32]. Over time, MDS have made immense contributions towards drug discoveries, exploring ligands such as small molecules, chemicals derived from plants, peptides, and proteins against targets such as protein kinases (PKs) [33][34][35][36][37], G-protein-coupled receptors (GPCRs) [38], and NMDA receptors [20,39,40]. Altogether, this method affords a means for drug design by providing a holistic approach to understanding the mechanism of receptor activation/deactivation, inhibiting the receptor to the mechanism of drug resistance. In the pharmaceutical industry, computational drug design has played a vital role in the discovery, design, and analysis of drugs. Computer technology nowadays is so rich and advanced that the accuracy of biomolecular simulations is consistently high enough to be used to truly drive preclinical drug discovery projects. Among the computational tools used for drug discovery, molecular dynamics simulations (MDS) and related methods are routinely used nowadays. Their main contribution is towards the understanding of structural flexibility together with entropic effects of complex systems. This allows an indepth estimation of the thermodynamics and kinetics related to drug-target recognition and binding [18]. MD simulations have been extensively used in the study of protein-protein and protein-ligand interactions, and the study of the mechanism of drug action [19][20][21][22][23]. The plethora of applications of molecular dynamics simulations extend from the study of complex and dynamic processes that play a central role in biological systems to the structure determination from X-ray crystallography and NMR experiments. In biological systems, the main application of MDS focuses on studies related to the stability of proteins through conformational changes and folding. MD studies also enable the molecular recognition of cellular components such as DNA, membranes along with complexes (drug-receptor), and ion transport [24][25][26][27]. MD studies have also enabled the study of the mechanism of drug resistance [28][29][30][31]. Molecular dynamics simulations calculate the binding free energy, which is very helpful in investigating receptor-ligand interactions [32]. Over time, MDS have made immense contributions towards drug discoveries, exploring ligands such as small molecules, chemicals derived from plants, peptides, and proteins against targets such as protein kinases (PKs) [33][34][35][36][37], G-protein-coupled receptors (GPCRs) [38], and NMDA receptors [20,39,40]. Altogether, this method affords a means for drug design by providing a holistic approach to understanding the mechanism of receptor activation/deactivation, inhibiting the receptor to the mechanism of drug resistance.
During MD simulations, configurations of the progressing system are sequentially generated; these configurations come out as trajectories that contain specific details of the positions and velocities of a particle over the time of simulation. These trajectories are explored to calculate a variety of properties, such as kinetic measures and free energy and other macroscopic quantities. These properties can be further compared with experimental observables that are helpful in drug design. This method was initially perceived within theoretical physics in the late 1950s, and its application was extended to chemical physics, materials science, biomolecular modelling, and, more recently, drug discovery [41].
In the past decade, there has been immense progress in the development of algorithms and technology using mutations, which has transformed the field of protein design and engineering to attain tailormade proteins suitable for pharmaceutical and biotechnological applications. Among various in silico tools available, the Cologne University Protein Stability Analysis Tool (CUPSAT) [42], Site Directed Mutator (SDM) [43], PopMusic 2.1 [44], SNPeffect 4.0 [45], PolyPhen-2 [46], DUET [47], MAESTRO Web [48], DynaMut [49], and mCSM PPI2 [50] for mutational studies have been successfully used to evaluate stability change (stabilizing or destabilizing) and, after mutations, to predict the phenotypic consequence of missense variants. These tools are structure based, sequence based, and energy based, and combined features (statistical approach and/or machine learning methods, such as neural networks and support vector machines (SVM). These methods are fast, user-friendly, and reliable and promise to be invaluable in the development of proteins with a wide range of impactful applications. The applications of these tools extend from understanding the origin of diseases caused by misregulation of protein maintenance [51,52] to discriminating disease-associated mutations from non-disease mutations, studying drugresistant mutations [53][54][55], and providing important structural and functional insights into designing new proteins [56][57][58][59]. To design DARPins as ERK2 inhibitors, a multi-point mutation approach, MAESTRO, was applied to the wild-type DARPin protein to identify the stabilizing mutation points, followed by validation of the binding energy of mutants employing MD simulations using MM-PBSA/GBSA protocols. The details of the results are discussed in the next section.

Design and Prediction of New Inhibitors
The effect of mutations on the thermodynamic stability of DARPins (E40) was analysed using MAESTRO and the predictions were further analysed using other algorithms, and a comparison of ∆∆G predictions is shown in Table S1. MAESTRO is an easy and standalone program that provides different kinds of mutation experiments on single chains and protein complexes. The predictive power of this method is suggested to be reliable as it combines multiple linear regression (MLR), a neural network approach (NN), and a support vector machine (SVM) that allows to include additional information such as protein size or solvent accessibility. The mutation sensitivity profile of E40 is shown in Figure 2. Mutation points were selected based on MAESTRO suggestions, from which only randomized residues of DARPin E40 were chosen to undergo further investigation; refer to Figure

Evaluation of Specific Mutations
In total, seven single-point favourable mutations (S380, I389, D421, N422, A443, D454, and R455) were suggested by MAESTRO based on predicted change in stability and confidence estimation calculated in terms of ∆∆G and Cpred, respectively. Wild-type amino acids (AA) in DARPin E40 were mutated with 16 other available amino acids (except cysteine, proline, and glycine). Table 1 shows the suggested mutations with their respective ∆∆Gpred and Cpred.

Evaluation of Specific Mutations
In total, seven single-point favourable mutations (S380, I389, D421, N422, A443, D454, and R455) were suggested by MAESTRO based on predicted change in stability and confidence estimation calculated in terms of ∆∆G and C pred , respectively. Wild-type amino acids (AA) in DARPin E40 were mutated with 16 other available amino acids (except cysteine, proline, and glycine). Table 1 shows the suggested mutations with their respective ∆∆G pred and C pred .

Evaluation of Selected Mutants by MDS
In the light of these results, out of all the suggested mutation points from MAESTRO, a total of 13 mutation points meeting the criterion, i.e., ∆∆G < 0.0 (stabilizing) and C pred values 1 (highly reliable), were selected for further analysis. Their structures were modelled using the "Build mutant" protocol in DS Modeller [60], further optimized, and finally subjected to large-scale MD simulations to investigate the structural consequences of mutating residues. To study the forces interplay, trajectories of mutant complexes were analysed for binding free energy, per-residue binding free energy of complexes, and pairwise binding free energy of residues within 4 Å.
The trajectories obtained from the production run of 100 ns MD simulations of mutants were analysed for their binding free energy. Binding free energies of mutant complexes (A443D/ERK2, S380L/ERK2, A443N/ERK2, and N422T/ERK2) suggested better binding than E40/ERK2 (Table 2. In the next section, the decomposition of binding energy and important interacting residues with their H-bonds are discussed in detail. MAESTRO gave the highest match with other prediction techniques and has been used as a guide for the new design. For ranking MM-GBSA, binding energy criteria are used and the energy contributions are G gas and G solv . The polar and non-polar contributions are EGB and ESURF, respectively, for the GB calculations shown in Table 3. Although the electrostatic contribution from all mutants is big, it is mostly compensated by a large positive polar contribution (EGB), making the total polar contribution (EEL + EGB) positive and hence it is mostly the van der Waals term that contributes towards the total binding free energy.
To understand the residue contribution from DARPin and ERK2, a comparison was made for the decomposed free energy for each residue. Important binding residues are shown with their respective energies for all the three mutants, and E40/ERK2.DARPin residues start with "L," while ERK2 residues start with "R" followed by a three-letter code of amino acids. From Figure 3 it is suggested that most of the important interactions from ERK2 come from the activation loop, specifically ARG180 contributing the most to the αG and L14 regions, while for all the DARPins (N3C), A443D, S380L, and D421W interactions come from all repeats (2,3), including L-ASP409, TRP413, and ASP421, but residues from C-cap (L-TYR444, ASP454, and PHE477) terminal contribute the most. Among all the complexes, A443D/ERK2 and S380L/ERK2 show stronger interactions with decomposed binding free energy <−5 kcal/mol for the above-mentioned residues.  Figure 3. Comparison between the per-residue free energy decomposition of E40/ERK2 and mutants. DARPin residues start with "L", while ERK2 residues start with "R" followed by a three-letter code of amino acids.

Exploring the Binding Mechanism of DARPins with ERK2
The main contribution towards the receptor-ligand affinity comes from non-covalent interactions. In the complex A443D/ERK2, alanine at position 443 on the DARPin loop is mutated to aspartate. Alanine contains a non-polar aliphatic R (CH3) group, which is small, being mutated to a polar aspartate that has a negatively charged or acidic R (CH2COO) group and a long side chain that can enhance electrostatic interactions, which can be seen from the electrostatic energy of −626.46 kcal/mol ( Table 3). The second-best selected mutant of E40 that came out after MDS is S380L with a −7 kcal/mol dip in binding free energy (−56.74 kcal/mol) compared to the previous one (Table 2). Serine contains a hydroxymethyl group and is classified as a polar amino acid mutated to a positively charged and non-polar leucine with a side chain containing an isobutyl group. This mutation leads to a more negative electrostatic interaction energy of −673.80 kcal/mol than E40/ERK2 (Table 3). Comparison between the per-residue free energy decomposition of E40/ERK2 and mutants. DARPin residues start with "L", while ERK2 residues start with "R" followed by a three-letter code of amino acids.
To evaluate the effect of stabilizing mutations on E40 DARPin, the RMSD values of the position differences of backbone atoms between mutant and wild-type structures were calculated throughout the MD simulations. A comparison of the RMSD of E40/ERK2 with the other four complexes pE59/ERK2, A443D/ERK2, S380L/ERK2, and D421W/ERK2 is shown in Figure S1. According to the figures, all the simulations are well converged. According to MAESTRO suggestion, the mutations A443D, S380L, and D421W are stabilizing. After performing 100 ns simulation on the complexes A443D/ERK2, S380L/ERK2, and D421W/ERK2, the average RMSD values are~3.19, 2.61, and 2.80 Å, respectively.

Exploring the Binding Mechanism of DARPins with ERK2
The main contribution towards the receptor-ligand affinity comes from non-covalent interactions. In the complex A443D/ERK2, alanine at position 443 on the DARPin loop is mutated to aspartate. Alanine contains a non-polar aliphatic R (CH 3 ) group, which is small, being mutated to a polar aspartate that has a negatively charged or acidic R (CH 2 COO) group and a long side chain that can enhance electrostatic interactions, which can be seen from the electrostatic energy of −626.46 kcal/mol ( Table 3). The second-best selected mutant of E40 that came out after MDS is S380L with a −7 kcal/mol dip in binding free energy (−56.74 kcal/mol) compared to the previous one (Table 2). Serine contains a hydroxymethyl group and is classified as a polar amino acid mutated to a positively charged and non-polar leucine with a side chain containing an isobutyl group. This mutation leads to a more negative electrostatic interaction energy of −673.80 kcal/mol than E40/ERK2 (Table 3).
To design a high binding inhibitor, it is required to understand interactions that differentiate it from the low binding inhibitors. For this purpose, the third mutant, D421W, having ∆∆G = −41.08 kcal/mol (Table 2), was selected. An aspartic acid at position 421 was mutated to tryptophan that contains a side chain indole, making it non-polar through its aromatic amino acid. Unlike the previous mutations that were located on the DARPin loops, this point mutation on the DARPin repeat has a negative effect on the binding affinity with a binding free energy of −41.08 kcal/mol compared to −49.50 kcal/mol for E40/ERK2. Figure 4 displays pairwise decomposition free energy for residues within 4 Å. The DARPin interacts with the receptor mainly through the activation loop, alpha G and MAPK insertion regions [61]. The receptor-ligand interactions are mainly through hydrogen bonds (details of the type of hydrogen bonds are presented in Table 4). For all the three mutants, the hotspot residues show interactions mostly through cation-Π and salt bridges. The two most remarkable salt bridge interactions occur between ARG180-ASP454 and LYS220-ASP409 ( Figure 5). In the mutant A443D, salt bridges ARG180-HH12-ASP454-OD1 and LYS220-HZ1-ASP409-OD2 possess quite a low polar interaction energy-(−8.76 and −8.25 kcal/mol) and (−8.98 and −8.27 kcal/mol), respectively; hence the total binding free energy of these residue pairs. There is a remarkable decrease in the energy of LYS220-HZ1-ASP409-OD2 here compared to E40/ERK2. For the next mutant complex S380L/ERK2, the pairwise interaction of 4 Å residues is almost like that in the previous case A443D/ERK2. Here, the interactions between receptor and ligand are also through hydrogen bonds (Table 5) but differ in their energies, with two strong salt bridge interactions between ARG180-HH22-ASP454-OD1 (−8.30 kcal/mol) and LYS220 HZ2-ASP409 OD2. A noticeable drop is observed in the decomposition energy of the receptor-ligand pair LYS220-ASP409 (−11.45 kcal/mol). For the third mutant complex, D421W/ERK2 also follows a pattern similar to that of the pairwise interaction of 4 Å residues like that in the previous cases of A443D and S380L. As in the previous cases, the interactions occur between receptor and ligand through hydrogen bonds, but are not as strong as the previous mutations. Except for a strong salt bridge interaction between ARG180-HH22-ASP454-OD1 (−8.30 kcal/mol), all other interactions have lower ∆∆G values compared to the other two mutants, while salt bridge LYS220 HZ2-ASP409 OD2 shows a remarkable increase (positive) in ∆∆G values. The main polar receptor-ligand interactions of the mutant complexes are shown in Figure 6. dues like that in the previous cases of A443D and S380L. As in the previous cases, the interactions occur between receptor and ligand through hydrogen bonds, but are not as strong as the previous mutations. Except for a strong salt bridge interaction between ARG180-HH22-ASP454-OD1 (−8.30 kcal/mol), all other interactions have lower ∆∆G values compared to the other two mutants, while salt bridge LYS220 HZ2-ASP409 OD2 shows a remarkable increase (positive) in ∆∆G values. The main polar receptor-ligand interactions of the mutant complexes are shown in Figure 6.        Electrostatic interactions are a pivotal player in protein interactions and helpful in understanding intermolecular protein-protein interactions as they are long-range and have influence on charge molecules. Electrostatic potential maps are also called electrostatic potential energy maps or molecular electrical potential surfaces. These maps demonstrate the charge distributions of molecules in three dimensions and allow us to visualize the variably charged regions of a molecule. Knowledge about the charge distributions can be useful in determining how molecules interact with one another. Moreover, electrostatic forces help in fast recognizing the right partner among hundreds of thousands of candidates present in the intracellular environment of protein-protein complexes [62].
To study in detail the specific electrostatic potential of the interacting residues, the interface region in the map was zoomed and the potential was matched with the van der Waals and electrostatic energies of the interacting residues in receptor-ligand pairs. In Figure 7, it can be observed that the interacting residues from receptors and ligands have opposite potentials, i.e., electropositive (blue) and electronegative (red). The electrostatic and van der Waals energy shown in the bottom left corner also corroborates the same. We know that van der Waals are weak forces or temporary attractions between electron-rich regions of one molecule and electron-poor regions of another. For the arginine residue from ARG180, a guanidino group is protonated to give the guanidinium form (-C-(NH 2 ) 2 + ), making arginine a charged, aliphatic amino acid; the lysine from LYS220 is also protonated −NH 3 + at physiological pH, and the tyrosine from TYR222 contains 4-hydroxyphenylalanine that is neutral. These residue charges extend from positive to neutral, which is shown by blue to white regions in the maps. In order to interact, they need to recognize oppositely (negatively) charged residues (aspartate form, −COO − ) in their vicinity that are ASP409, ASP454, and ASP475, shown by red regions in the maps. ARG180-ASP454 shows the strongest electrostatic energy because of the salt bridge formed In general, single-point mutation in ankyrin does not affect the secondary structure. We have predicted the protein secondary structure for the wild-type and the mutated ankyrin using PredictProtein (ROSTLAB, Technische Universität München, https://predictprotein.org/) [64] and SCRATCH protein predictor (Institute for Genomics and Bioinformatics, University of California, Irvine, USA, http://scratch.proteomics.ics.uci.edu/) [65]. Results are shown in Supplementary Material, Figure S2. It was observed that the secondary structures of mutant DARPins were unaffected.

Discussion
To design DARPins with higher binding affinity towards ERK2, the MAESTRO method was used to predict stability upon point mutations in N3C DARPin (E40). From the suggested mutation points, 13 stabilizing mutants were subjected to 100 ns simulations and then compared with the wild-type complex (E40/ERK2) in terms of binding free energy. After evaluation, two DARPins, A443D and S380L (showing higher binding affinity than E40), and D421W (showing lower binding than E40) were selected to further analyse their binding mechanism. Out of the suggested mutants, the key elements of the ankyrin mutants with high and low affinity towards ERK2 were compared in terms of binding free energy, decomposition free energy, and strength and type of hydrogen bonds formed between receptor and ligands. It was found that all the DARPins show interactions, mostly coming from their third repeat and C-cap terminal residues. The most persistent interactions in each structure were studied via a H-bond analysis, wherein it was found that the high binding affinity of all the N3C DARPins is mainly attributed to salt bridges (ARG180-ASP454 and LYS220-ASP409) along with various other DARPin-ERK2 interactions. Analysis of the energy decomposition shows that electrostatic and van der Waals interactions were the main contributors towards the total binding free energy of In general, single-point mutation in ankyrin does not affect the secondary structure. We have predicted the protein secondary structure for the wild-type and the mutated ankyrin using PredictProtein (ROSTLAB, Technische Universität München, https: //predictprotein.org/) [64] and SCRATCH protein predictor (Institute for Genomics and Bioinformatics, University of California, Irvine, USA, http://scratch.proteomics.ics.uci. edu/) [65]. Results are shown in Supplementary Material, Figure S2. It was observed that the secondary structures of mutant DARPins were unaffected.

Discussion
To design DARPins with higher binding affinity towards ERK2, the MAESTRO method was used to predict stability upon point mutations in N3C DARPin (E40). From the suggested mutation points, 13 stabilizing mutants were subjected to 100 ns simulations and then compared with the wild-type complex (E40/ERK2) in terms of binding free energy. After evaluation, two DARPins, A443D and S380L (showing higher binding affinity than E40), and D421W (showing lower binding than E40) were selected to further analyse their binding mechanism. Out of the suggested mutants, the key elements of the ankyrin mutants with high and low affinity towards ERK2 were compared in terms of binding free energy, decomposition free energy, and strength and type of hydrogen bonds formed between receptor and ligands. It was found that all the DARPins show interactions, mostly coming from their third repeat and C-cap terminal residues. The most persistent interactions in each structure were studied via a H-bond analysis, wherein it was found that the high binding affinity of all the N3C DARPins is mainly attributed to salt bridges (ARG180-ASP454 and LYS220-ASP409) along with various other DARPin-ERK2 interactions. Analysis of the energy decomposition shows that electrostatic and van der Waals interactions were the main contributors towards the total binding free energy of these systems. The binding free energy of the newly designed DARPins was improved up to~20%. The binding affinity of these systems shows a trend of E40 < D421W < S380L < A443D.
The wild-type DARPin (E40) along with mutants interact (A443D, S380L, and D421W) through hydrogen bonds, mainly salt bridges and cation-Π. Moreover, the strength of the salt bridges formed between ARG180-ASP454 and LYS220-ASP409 plays a significant role in deciding the binding affinity of DARPins towards ERK2. For a DARPin to have a good binding affinity, both interactions must be strong, which can be shown by the order of decomposition free energy of these receptor-ligand pairs and the total binding free energy of the complex. Electrostatic potential surface analysis revealed that mutants that can generate an electronegative potential near the binding interface, show good binding with ERK2. It was also observed that mutations on the loops of DARPins extend better binding compared to that on a DARPin repeat. All this information can be used in the design of new DARPin inhibitors against ERK2.

Materials and Methods
The starting structure for the present study is the X-ray crystal structure of DARPin E40 complexed with ERK2 (PDB ID: 3ZU7) [66], designated as E40/ERK2 for further reference in this study. For the complex (E40/ERK2), chain A of ERK2 and chain C of DARPin (E40) were taken as the initial structure for the MD simulations ( Figure 8). these systems. The binding free energy of the newly designed DARPins was improved up to ~20%. The binding affinity of these systems shows a trend of E40 < D421W < S380L < A443D.
The wild-type DARPin (E40) along with mutants interact (A443D, S380L, and D421W) through hydrogen bonds, mainly salt bridges and cation-Π. Moreover, the strength of the salt bridges formed between ARG180-ASP454 and LYS220-ASP409 plays a significant role in deciding the binding affinity of DARPins towards ERK2. For a DAR-Pin to have a good binding affinity, both interactions must be strong, which can be shown by the order of decomposition free energy of these receptor-ligand pairs and the total binding free energy of the complex. Electrostatic potential surface analysis revealed that mutants that can generate an electronegative potential near the binding interface, show good binding with ERK2. It was also observed that mutations on the loops of DARPins extend better binding compared to that on a DARPin repeat. All this information can be used in the design of new DARPin inhibitors against ERK2.

Materials and Methods
The starting structure for the present study is the X-ray crystal structure of DARPin E40 complexed with ERK2 (PDB ID: 3ZU7) [66], designated as E40/ERK2 for further reference in this study. For the complex (E40/ERK2), chain A of ERK2 and chain C of DARPin (E40) were taken as the initial structure for the MD simulations ( Figure 8). MD simulations were performed using the GPU version of Particle Mesh Ewald Molecular Dynamics (PMEMD.CUDA) from AMBER14 [67]. At the molecular level, physical forces were implemented using ff14SB [68] protein force field to carry out MD simulations. The structures were solvated (using the tLeap module implemented in AMBER). The water model TIP3P [69] was used, wherein a cubic box of water extends at least 10 Å from the solute in each direction. A cut-off distance of 15 Å was used to compute the nonbonded interactions (electrostatic interactions and van der Waals interactions). For each of the complexes, a 100 ns long simulation was performed. To minimize the edge effects, all simulations were performed under periodic boundary conditions, and to treat longrange electrostatics, the particle mesh Ewald method was used [70]. To relax the system prior to MD simulation, the complexes were minimized using a series of steepest descent (SD) and conjugated gradient (CG) under the sander module of the AMBER14 program. During the simulation, the system was heated gradually over a period of 60 ps from 0 to 310 K (biological temperature), and a force constant of 5 kcal/mol Å 2 was applied to restrain the atomic position. A constant pressure of 1 atm for 200 ps (NPT) was applied MD simulations were performed using the GPU version of Particle Mesh Ewald Molecular Dynamics (PMEMD.CUDA) from AMBER14 [67]. At the molecular level, physical forces were implemented using ff14SB [68] protein force field to carry out MD simulations. The structures were solvated (using the tLeap module implemented in AMBER). The water model TIP3P [69] was used, wherein a cubic box of water extends at least 10 Å from the solute in each direction. A cut-off distance of 15 Å was used to compute the non-bonded interactions (electrostatic interactions and van der Waals interactions). For each of the complexes, a 100 ns long simulation was performed. To minimize the edge effects, all simulations were performed under periodic boundary conditions, and to treat long-range electrostatics, the particle mesh Ewald method was used [70]. To relax the system prior to MD simulation, the complexes were minimized using a series of steepest descent (SD) and conjugated gradient (CG) under the sander module of the AMBER14 program. During the simulation, the system was heated gradually over a period of 60 ps from 0 to 310 K (biological temperature), and a force constant of 5 kcal/mol Å 2 was applied to restrain the atomic position. A constant pressure of 1 atm for 200 ps (NPT) was applied under Langevin dynamics, followed by a 40 ps volume-constant period (NVT) at a force constant of 2.5 kcal/mol Å 2 , which was maintained and followed by 100 ps dynamics at a force constant of 1.25 kcal/mol Å 2 . Finally, unrestrained production runs were performed for 100 ns, wherein no force was applied on any protein atoms in the NVT ensemble at a constant temperature of 310 K (biological temperature). For all analyses, 500 snapshots were taken from the last 5 ns of the simulation (96-100 ns). To check the system equilibrium, root-mean-square deviation (RMSD) of all backbone atoms was performed using the cpptraj module [71] incorporated in AmberTools 15 and compared to the starting structure. All simulations were carried out under periodic boundary conditions [72]. To treat long-range electrostatics, the particle mesh Ewald method was used [73][74][75], and the SHAKE algorithm was employed to constrain bond lengths involving hydrogen atoms. A 2 fs time step was set up while the trajectory was recorded every 0.1 ps. To relax the system prior to MD simulation, a series of steepest descent (SD) along with conjugated gradient (CG) minimizations with a total of 500 steps each was performed. Post-processing of the trajectories was performed using the MM-PBSA/GBSA protocol [76].
Supplementary Materials: The following are available online: Figure S1: Comparison of the RMSD values of mutants and E40/ERK2; Figure S2: Comparison of secondary structures of DARPin E40 and mutants (A443D, S380L, and D421W) obtained by using (A) PredictProtein and (B) SCRATCH protein predictor web servers. (For PredictProtein, helixes are shown in blue and others in yellow; for the SCRATCH protein predictor, H is for helix and C is for others.); Table S1: Comparison of predicted ∆∆G (kcal/mol) mutants of 3ZU7 from different web-based algorithms.

Data Availability Statement:
The data presented in this study are available within the article or supplementary materials.