Various NMR buffer conditions were explored, to determine conditions which gave maximum stability for NMR studies, and highest quality 2D HSQC spectra. The conditions examined were deuterated Tris buffer and potassium phosphate buffer, with pH values of 5.0 and 6.8. The procedure outlined below is for the 50 mM Tris, 100 mM KCl, and 2 mM DTT buffer samples at pH 5.0 and 6.8. This same procedure was followed for the investigation in buffer containing 20 mM potassium phosphate, 75 mM KCl, and 2 mM DTT.

An ¹⁵N labeled DUSP5pdRDB sample with a concentration of 500 µM was thawed on ice from −80°C storage. Samples without glycerol-d₈ were prepared by transferring 300 μL of the ¹⁵N labeled protein and 30 μL D₂O to a 1.5 mL Eppendorf tube. Glycerol-d₈ containing samples were prepared in a similar fashion by transferring 300 μL of ¹⁵N labeled protein, 30 μL D₂O, and 30 μL glycerol-d₈ for a final concentration of 10% D₂O and 10% glycerol-d₈. The samples were spun down in a centrifuge at 14,000 rpm for 2 min to remove any particulate. The supernatant was transferred to a Shigemi NMR tube.

Hydrolysis of the pThr-Glu-pTyr peptide was monitored using ³¹P 1D NMR and a Varian 500 MHz spectrometer equipped with a broad band probe. 1 mM of peptide was incubated with 100 μM of DUSPpdWT protein for 3 min, and spectra were collected at 25^oC. See Figure 2 and, for a time course, Supplementary Figure S3.

2D ¹H-¹⁵N HSQC spectra were collected using a 500 MHz Varian NMR spectrometer of a sample of expressed and purified DUSP5pdRDB protein labeled with ¹⁵N and ¹³C. The temperature for the experiments was adjusted to either 25°C, 30°C, or 37°C. After the sample was done running on the NMR spectrometer it was transferred to an oven set to either 25°C, 30°C, or 37°C. The samples were held at their respective temperatures for a total of 72 h 2D ¹H-¹⁵N HSQC experiments were run on each sample after the completion of each 24 h period until 72 h. Final optimized conditions, based on quality of 2D ¹H-¹⁵N HSQC spectra, were determined to be 20 mM potassium phosphate (pH 5), 75 mM KCl, 2 mM DTT with 10% glycerol, with NMR experiments run at 25^oC. These are the conditions used for the 3D NMR experiments, and for chemical assignments, using the DUSP5pdRGB construct (Supplementary Figure S1).

Using these optimized conditions, protein chemical shift assignments were made for ¹⁵N and ¹³C labeled and cleaved DUSP5pdRDB protein (summarized in Supplementary Table S1). The multidimensional NMR experiments were performed on a 600 MHz Bruker NMR spectrometer equipped with a cryoprobe, using instrumentation at Medical College of Wisconsin.

The sample was concentrated to 1 mM in a total volume of 0.3 mL. 30µL glycerol-d₈ was added to the sample and it was flash frozen and stored in the −80°C freezer until ready to analyze with the NMR spectrometer. 30µL D₂O and 0.02% sodium azide were added to the sample. The sample was transferred to a Shigemi NMR tube and a standard suite of 2D and 3D experiments were collected for structure determination: 2D ¹H-¹⁵N HSQC, 3D HNCO, HN(CO)CA, HN(CO)CACB, CCONH, HN(CA)CO, HNCA, and HNCACB.

Spectral data were processed using NMRpipe, then experiments were viewed in Xeasy and were peak picked. Once all the peaks were picked in all of the experiments, the data were run through Garant, which was able to auto assign roughly 85% of the chemical shifts. Next, the Garant assignments were opened into the HSQC experiment in Xeasy. The assignments were checked and corrected if necessary. This procedure was repeated for the HNCO experiment. The correct peak list of the HNCO experiment were loaded into the HN(CO)CA. The corrected HN(CO)CA list was loaded into the HNCACB, and the HN(CO)CACB experiments. Lastly, the HN(CO)CACB peak list was loaded into the CCONH experiment. The completed peak list was once again opened into the HSQC file for a final check of chemical shifts. These chemical shifts were saved to the prot file, with final assignments provided in Supplementary Table S1. Assignments mapped onto 2D ¹H-¹⁵N HSQC spectra are in Supplementary Figure S4, S5.

The chemical shifts that were assigned for the primary structure were re-formatted for the TALOS+ web based secondary structure assignment (National Institute of Diabetes and Digestive and Kidney Diseases, 2014). Once the file was reformatted, the data file was uploaded to the TALOS+ server (Shen et al., 2009). This produced the analysis files. These files were then read through the jRAMA+ web-based java viewer for TALOS+ torsion angle predictions. These predictions were compared to the crystal structure, with pdb code 2G6Z (See Figure 3).

After dialysis, two NMR sample tubes were prepared by adding 275 µL of 564 µM DUSP5pdRDB protein, 30 µL glycerol-d₈, 30 µL D₂O, and varying concentrations of the Ac-pThr-Glu-pTyr-NH₂ (pT-E-pY) peptide (N-terminus acetylated and C-terminus amide) derived from pERK. Peptides were synthesized (Biomatik LLC, Wilmington, Delaware) to greater than 98% purity as determined by HPLC. The concentration of the tripeptide added was 0.250 mM. The Eppendorf tubes were centrifuged at 14,000 rpm for 2 min. After centrifugation the supernatant of each Eppendorf was transferred to each of two Shigemi NMR tubes. 2D ¹H-¹⁵N HSQC spectra were collected on a 500 MHz Varian NMR spectrometer, at 25°C. Spectra were processed then overlaid using vNMRDraw software (shown in Figure 4). Exchange-broadened active site residues (identified in Figure 4) as well as residues with low predicted order parameters (low S²; from Figure 3) are indicated in our previously reported model of the DUSP5-ERK2 complex (Figure 5).

Peptide stock solutions were kept in buffer. Steady state kinetics were performed as described for the peptide inhibition studies, but monitoring actual dephosphorylation of the three tripeptides by DUSP5pdWT: pT-E-Y, T-E-pY, pT-E-pY. 1 mM of each tripeptide was incubated with 0.18 mM DUSP5pdWT enzyme at 25^oC, and at timed intervals samples were analyzed for inorganic phosphate using the Biomol Green reagent (at 620 nm) (Supplementary Figure S6). Lack of protein aggregation at these concentrations is evidenced by consistently narrow linewidths in 2D HSQC spectra at various enzyme concentrations (i.e., minimal T2 broadening in Figs. 4, S4 and S5) as well as our previous data on dynamic light scattering of DUSP5 (Gupta et al., 2019). Kinetic studies were also repeated with varied tripeptide substrate concentration, to generate Michaelis-Menten curves. Results of fitting to Equation (1) are shown in Table 1. v = V max S K m + S (1)

Recombinant DUSP5 PD (DUSP5pdWT) phosphatase activity was assayed as previously described, with minor modifications (Nayak et al., 2014; Neumann et al., 2015). All assays were performed in clear bottom 96-well microplates (Greiner 655001) with a total assay volume of 210 μL, using a SpectraMax M5 microplate reader (Molecular Devices). Assay buffer contained 100 mM Tris, 100 mM NaCl, 5 mM MgCl₂·6H₂O with 1 mM dithiothreitol (DTT) at pH 7.5. p-nitrophenol phosphate (pNPP, Sigma Aldrich) at 5 mM, a concentration near the established K _m , was used for substrate. DUSP5 PD-mediated hydrolysis of pNPP to p-nitrophenolate, a chromogenic product which absorbs at 405 nm with a molar extinction coefficient of 18,000 M⁻¹ cm⁻¹, was monitored over time using the Spectramax M5 microplate reader. Wells that contained assay buffer with pNPP were used as blanks. Stock solutions of peptides were kept in buffer. Serial dilutions of each stock solution were prepared to generate a range of inhibitor concentrations for the steady state inhibition studies. Absorbance values were converted to p-nitrophenolate concentration using the product extinction coefficient and a micro well path length of 0.545 cm DUSP5 PD (DUSP5pdWT) enzyme progress curves over a range of substrate and peptide inhibitor concentrations were fitted globally using GraphPad Prism 6 to Equation (2) for steady-state competitive inhibition: v = V max S K m 1 + I K i + S (2) where v is the initial velocity, V max is the maximum velocity, K m is the Michaelis constant, S is the pNPP concentration (varied from 1 to 80 mM), I is the peptide inhibitor concentration, and K i the fitted inhibition constant. The parameters V max , K m , and K i are shared in a global fit, resulting in one best-fit value for the complete data set (Table 2; Supplementary Figure S7; Supplementary Table S2).

All assays were carried out using pNPP as a substrate in the presence of DTT using a microplate reader and a total volume of 100 μL. Each assay contained an enzyme concentration of approximately 0.072 mg/mL (appropriately 3.13 μM DUSP5pdWT) which was left in a stock solution at about 1 h at 4°C before use. Each assay consisted of a reaction mixture in 100 mM tris-hydrochloride buffer containing 100 mM NaCl and 5.0 mM magnesium chloride (MgCl₂) at pH 7.5. The concentration of pNPP was varied from 2 to 25 mM (2, 4, 8, 15 and 25 mM) which was determined after measuring a K_m value of 8.02 mM using these conditions. The kinetic parameters were measured in the presence of 0.045 mM DTT in the assay solution and the assay was started with the addition of pNPP.

The inhibition by vanadate was carried out using inhibitor concentrations of 100, 200, 400, 600 and 800 nM under the assay conditions described above and in an assay solution containing 0.045 mM DTT. For phosphate inhibition studies, the phosphate inhibitor concentrations used were 0.5; 1; 2; 5; 10 mM. Both sets of inhibition studies are shown in Supplementary Figure S8, with both phosphate and vanadate behaving as competitive inhibitors. The kinetic runs were repeated three times and the data were initial processed using excel and the global kinetics analysis was performed using GraphPad Prism 10.1.2. The errors provided in the analysis were calculated by the program.

Three series of molecular dynamics simulations with varied structural features of the molecules were conducted: neutral cysteine (Cys263) in the active site of DUSP5 and phosphorylated residues of the tripeptides in dianionic PO₃ ^2- form (Series 1), deprotonated Cys263 and phosphorylated residues of the tripeptides in dianionic PO₃ ^2- form (Series 2), deprotonated Cys263 and phosphorylated residues of the tripeptides in monoanionic PO₃H⁻ form (Series 3). Results are summarized in Table 4 and Supplementary Figure S9, S10.

Calculations were done using methods described previously (Talipov et al., 2016; Gupta et al., 2019). Briefly, MD simulations of DUSP5pdWT were performed with the secondary binding site and active/catalytic binding sites loaded or unloaded with mono-phosphorylated ERK2 (i.e., pT-ERK2 or pY-ERK2) tripeptides described above (T-E-pY; pT-E-Y; also T-E-Y), while the remaining site was vacant and conformationally flexible in both cases. These simulations of the PD with the loaded secondary binding site or active site serve to probe the structural rearrangement of the active and secondary binding sites as compared with the original X-ray structure of PD domain. MD simulations were performed as before (Gupta et al., 2019) using the Amber MD software (Case et al., 2005) with the ff14SB force field with a 10.0 Å force cutoff and the Particle Mesh Ewald algorithm (Duan et al., 2003) to treat long-range electrostatic interactions. A triplicate 8 microsecond MD simulation (24 μs in total) was performed with each of the tripeptides: pT-E-pY, pT-E-Y, T-E-pY, and T-E-Y. Solvent was modeled explicitly with the TIP3P water model (Jorgensen et al., 1983). For each simulation, the initial structure (prepared as described below) was subjected to 500 steepest descent (SD) + 500 conjugate gradient (CG) minimization, in which protein structure was frozen, and subsequent 1000 SD + 1500 CG full minimization steps performed. The resulting structure was subjected to a warm-up NVT simulation with temperature being increased from 0 to 300K during 40 ps and subsequent 3-ns NPT equilibration. The sets of coordinates and velocities after each nanosecond of simulation were used as a starting point for three independent production run simulations. Production runs were performed using an NPT ensemble, with the time constant for heat bath coupling for the system equal to 10 ps and a Monte Carlo barostat. For each production run simulation, it was confirmed that temperature, mass density, and potential energy were nearly constant and did not show any noticeable drift.

A docking study with molecular dynamics, as described above, was performed also using the fully phosphorylated peptide (pT-E-pY), and a population analysis was used to assess which of the two phosphate binding pockets was preferentially occupied by pT or pY, and at what energy (Figure 6).

The starting structure for simulations was the 2.7 Å resolution X-ray structure (PDB code 2G6Z) (Jeong et al., 2006), with Ser263 typically replaced with Cys, which is the WT (catalytically active) form. Structures were prepared using UCSF Chimera (Pettersen et al., 2004). Statistical analysis of the conformational dynamics was performed using R software version 3.3.1. For the cluster analysis, a hierarchical cluster analysis technique was used for the combined set of MD trajectories.

DUSP5 is a dual specificity phosphatase that cleaves both phospho-threonine (pT) and phospho-tyrosine (pY) residues, and it does this via a two-pocket gated mechanism that we have previously described (Figure 1). The peptide from the ERK activation loop that is dephosphorylated by DUSP5 actually contains proximal phosphates on both of these residues, separated by a glutamate: pT-E-pY. Studies presented herein characterize the binding and reaction of this tripeptide, derived from ERK. This pT-E-pY tripeptide is a substrate for DUSPpdWT with both phosphates cleaved 100%, based on ³¹P NMR analysis of the reaction progress (Figure 2). This result demonstrates that both the phosphate on the pY and the pT residues can be hydrolyzed by DUSPpdWT.

Hydrolysis was also examined at 10-fold lower enzyme concentration (Supplementary Figure S3) with the objective of observing a preference for either the pT or the pY residue. At the first timepoint after 3 min incubation, about 50% conversion was observed. Thus, about 25% each of the pT and pY peptides remained even after prolonged incubation. While both phosphates (pY and pT) are hydrolyzed efficiently in ³¹P NMR studies at 10-fold higher enzyme concentration, the fact that no more hydrolysis is occurring at lower enzyme concentration suggests either the enzyme is decomposed, or the enzyme is inhibited. In addition, the question remains as to whether there is differential affinity or binding order for pT or pY in each of the two binding pockets of DUSP5 (Figure 1), and how the enzyme catalyzed reaction occurs under steady state conditions.

Based on DUSP5 chemical shift values, secondary structure was predicted, with these predictions mapped onto the known DUSP5 x-ray crystal structure in Figure 3B. NMR data confirm that the solution structure is the same or at least consistent with the crystal structure, with perhaps more flexibility (loop-like structure) observed in some regions at the end of helices (Figure 3). Furthermore, the NMR data yielded a prediction of high flexibility (low S² values) in regions 215-220 and 232-237, which we term “Flexible Regions A and B.” Extensive exchange broadening in 2D ¹H-¹⁵N HSQC spectra (Figure 4) is consistent with protein conformational mobility in general. Mapping the Flexible Regions onto our modeled DUSP5-ERK complex reveals that these two loop regions define the binding interface with ERK (Figure 5), perhaps positioning the DUSP5 phosphatase domain relative to ERK so that the ERK activation loop (containing the pT-E-pY peptide) fits into the DUSP5 active site pocket for hydrolysis to occur.

Chemical shift perturbation studies, comparing 2D ¹H-¹⁵N HSQC spectra before and after binding of the pT-E-pY tripeptide, indicate there is significant exchange broadening of cross peaks (Figure 4), including primary active site residues Ala265, Ile267 and Ser268. This suggests that the DUSP5-(pT-E-pY) complex is in a dynamic rather than rigid state, consistent with our previously published mechanism (Figure 1). These residues are shown on the DUSP5 phosphatase domain structure in Figure 6A and are close in space to the catalytic Cys163 residue in the primary active site pocket. Their flexibility on the NMR exchange timescale (typically msec) is consistent with motions and flexibility in the primary active site pocket that may be necessary to transfer pY or pT from binding in one pocket to the next.

Steady state kinetic studies were performed with the tripeptides acting as substrates (vs inhibitors), with progress curves and Michaelis-Menten plots shown in Supplementary Figure S6. Fitted Michaelis-Menten kinetic profiles yielded a relative activity as a substrate, expressed as k_cat/K_m, of: pT-E-pY > T-E-pY > pT-E-Y (Table 1). This clear trend is only with the kinetic parameter k_cat/K_m, the “catalytic efficiency,” that reflects all steps from (and including) initial binding up to the first irreversible step (Figure 1C), which is most likely the (initial) phosphate transfer that creates the phosphor-enzyme adduct (as also reported by Brandão, et al. (2010), for PTP1B). What is not clear from these results is whether the pT-E-pY peptide binds sequentially in one site then the next, and—even though there is higher overall affinity for pY over pT—what the preference is in each of the two binding pockets. In other words, the binding mechanism needs further exploration, which we performed below using kinetic studies with pNPP as substrate and phosphorylated peptides and with molecular modeling.

Steady state inhibition studies were performed with pNPP as substrate and the fully phosphorylated, nonphosphorylated, and both mono-phosphorylated T-E-Y tripeptides as inhibitors. The absorbances at 405 nm, measuring the dephosphorylation of pNPP, was analyzed with global fit analysis fitting to the equation for competitive inhibition (Supplementary Figure S7). The resulting kinetic parameters (K_m, k_cat, K_i and k_cat/K_m) summarized in Table 2 indicate an overall binding affinity preference and relative phospho-tripeptide inhibition potency of: pT-E-pY > T-E-pY > pT-E-Y > T-E-Y. The diphosphorylated peptide is a roughly 3x more potent binder than the pY monophosphoryled peptide, which is in turn 3x more potent than the pT monophosphorylated peptide. But the greatest affinity boost comes from the addition of the initial phosphate group, since the pT monophosphoryleted peptide binds 30x stronger than nonphosphorylated peptide.

These data demonstrate that the phosphorylated peptides are inhibitory. While both pNPP and the phosphorylated peptide inhibitor can hydrolyze during the course of the assay, the rate of hydrolysis of pNPP is much faster than of the phosphorylated peptide. This can be rationalized because the pNPP concentration range is 2–20 mM whereas the peptides are present at much less than 1 mM, and the amount of the latter hydrolyzed is not contributing significantly to the overall level and rate of phosphate hydrolysis. Hence incubation of phosphorylated peptides with DUSP5pdWT will inhibit this enzyme and the observed ³¹P NMR results with low enzyme concentration where the reaction stops at 50% conversion can be explained at least in part by incubation of the phosphorylated peptide. Furthermore, these data show that the pY binds tighter than pT, consistent with the possibility that the initial binding of the diphosphorylated peptide will be through binding of pY. This hypothesis is explored further in molecule docking and dynamics experiments described below.

Since DUSP5 has an active site cysteine that is thought to act as a nucleophile (Figure 1), forming a Cys-phospho intermediate analogous to PTP1B, it might be expected DUSP5 is potently inhibited by vanadate, since Brandão, et al. (2009 and Brandao et al., 2010) had previously reported such a potent inhibitory adduct formed between vanadate and the cysteine thiol of PTP1B. Vanadate readily forms trigonal bipyramidal structures and are widely accepted as mimicking transition state or high energy intermediate structures for phosphate ester hydrolysis (Crans et al., 2014; Crans, 2015; McLauchlan et al., 2015). The structural advantage that vanadate has over other analogs including tungstate (WO₄ ^2-), molybdate (MoO₄ ^2-), magnesium fluoride (MgF₃ ⁻), and nitrate (NO₃ ⁻) is that it has the same overall charge as phosphate. Herein we explored the inhibition by vanadate and compared it with inhibition of phosphate, as our lab has done previously (Crans, 2015; McLauchlan et al., 2015). However, because of the redox properties of vanadate (Crans et al., 2004), it can be challenging to use in the presence of a reducing agent (DTT) that is needed to maintain the stability of protein phosphatase, with a readily oxidized cysteine in the active site (Crans et al., 2014; Crans, 2015; McLauchlan et al., 2015). We carried out preliminary studies with DUSP5 using varying concentrations of DTT. Tris forms a weak complex with vanadate, and therefore was compared to Hepes, which is a buffer with similar pK_a value that interacts minimally with vanadate (Crans et al., 1989). Preliminary studies in both buffers showed that some variations in K_m values for pNPP and K_i values for vanadate were observed. These preliminary studies allowed us to choose Tris buffer with minimum DTT concentration for studies with the DUSP5pdWT enzyme for comparison with the kinetic analysis carried out in this work.

Using the conditions determined by the pNPP K_m of 8.02 mM that was measured and calculated using global kinetics analysis, with vanadate at concentrations at 100, 200, 400, 600 and 800 nM, a K_i value of 638 nM was obtained for vanadate (Supplementary Figure S8). Considering that the active site of DUSP5 contains a cysteine residue and it has a secondary binding site which does not, DUSP5 is distinct from than other cysteine containing phosphatases. However, as shown by our data, we find that the inhibition is competitive and the K_i value is very potent, submicromolar. A summary of K_i values surveyed from the literature is listed in Table 3, and demonstrates the variations in K_i values reported for vanadate, and for comparison phosphate, considering the differences in phosphatases with cysteine in the active site, substrate, and pH values in the assay (Lawrence and Van Etten, 1981; Zhou et al., 1993; Huyer et al., 1997; Zhang et al., 1997; Zhang. and Zhang, 1998; Deng et al., 2002; Heo et al., 2002; Lu et al., 2012). Except for one bovine liver enzyme (Lawrence and Van Etten, 1981) these phosphatases all have micromolar or submicromolar K_i values. The most well studied vanadate-phosphatase system with a cysteine residue in the active site is protein phosphatase 1B (PTP1B) (Brandão et al., 2009; Brandão, et al., 2010; Crans, 2015). There are multiple X-ray structures of the vanadate-PTP1B complex showing a range of different complexes. One complex in particular has part of the peptide tyrosine-phosphate substrate bound to the PTP1B. The substrate is a hexapeptide with a sequence DADEYL which is a fragment of the epidermal growth factor receptor, which, in phosphorylated form (DADEpYL) exhibits strong affinity for PTP1B (Pedersen et al., 2004). A previous crystallographic study shows a Michaelis complex between a catalytically inactive C215S mutant of PTP1B and the phosphorylated peptide (Jia et al., 1995) but the second transition state of this protein with vanadate (Pedersen et al., 2004; Brandão, et al., 2010) has lost the tyrosine peptide. The crystal structure of the C215S mutant of PTP1B and the phosphorylated peptide shows the peptide extending from the active site to the adjacent site which provides a number additional stabilizing H-bonds and as such a system to which there is more direct comparison to DUSP5 than the other Cys-containing phosphatases.

For comparison, the inhibitory effect of phosphate on DUSP5 was also determined (Supplementary Figure S8B). Phosphate also exhibits a competitive inhibition pattern with DUSP5 and yields a K_i of 17.3 mM, similar to those reported previously (Lawrence and Van Etten, 1981; Zhang et al., 1997; Zhang. and Zhang, 1998). The considerably lower inhibitory effects of phosphate (5 orders of magnitude higher K_i) compared to vanadate (K_i = 638 nM) is due to the fact that phosphate is not able to form the stable pentacoordinate complex in the trigonal bipyramidal structure that vanadate can adopt, resembling the transition state or high energy intermediate as reported for other protein phosphatases (McLauchlan et al., 2015). Our studies demonstrate that DUSP5 is subjected to the same effects during its catalytic cycle, and in this regard resemble PTP1B.

Binding of the pT-E-pY tripeptide, and its monophosphorylated analogs, was explored using docking and molecular dynamics to assess preference for the two DUSP5 phosphate binding pockets (catalytic and the secondary site; Figure 1) and to verify if the pY residue binds tighter than the pT residue. A 24-microsecond molecular dynamics (MD) simulation of the DUSP5 complex with the pT-E-pY tripeptide was performed, and which showed that (a) the complex is stable and maintains pY in the catalytic site and pT in the secondary site (Figure 6A), (b) there is preferential binding of the pY group of pT-E-pY before pT (Figure 6B), (c) distance between the active site residue (Cys263) is shorter to pY than to pT (Figure 6C), and (d) the pY and pT groups bind in the catalytic (active) and the secondary binding site, respectively, with binding of pT in the secondary pocket permitting binding of pY in the active site, consistent with our published mechanism (Figure 1; Talipov et al., 2016). The results from the accelerated sampling of free energy profile for the complex formation (replica exchange umbrella sampling technique) are also consistent with the catalytic efficiency measurements: pT-E-pY (∆G = −6.8 kcal/mol) > T-E-pY (−6.5 kcal/mol) > pT-E-Y (−6.3 kcal/mol) > T-E-Y (−5.9 kcal/mol).

Our calculations show that the active site (containing the catalytic residue Cys263) is preferred over the secondary active site for binding phosphorylated residues. Tripeptides pT-E-Y and T-E-pY bind to the active site (∼5 Å) by their single phosphate groups in more than 50% of frames in all series in the MD simulation. pY is more likely to be involved in this interaction than pT. Numbers of frames with corresponding interactions for T-E-pY and pT-E-Y are about 80–90% and 50–60% respectively. Among all the series of simulations, the strongest interactions are observed in Series 3 with deprotonated Cys263 and the phosphate groups in monoanionic PO₃H⁻ form; and, Series 2 with deprotonated Cys263 and phosphate groups in dianionic PO₃ ^2- form showed poorest binding between modeled fragments, due to high negative charge of the ligands (Table 4; Supplementary Figure S9, S10).

In the case of a tripeptide with two phosphate groups (pT-E-pY), the main active site is occupied by the pY residue most of the time. The most illustrious results in this connection are from Series 3: the distribution between frames with Cys-pY and Cys-pT interactions are about 90% and 15% respectively (in 5% of cases both phosphorylated groups are located near Cys). According to the Boltzmann distribution, this ratio between conformers corresponds to 1.3 kcal/mol of difference in free energy. Calculations from Series 2 show that two dianions make the ligand too negative for interaction with Cys^-. In this case the active site in most of the frames is empty, while the occupation of the secondary binding site is considerably higher than in other sets of calculations. In simulations involving neutral cysteine, pT did not show interactions with the main active site, while pY is bound proximal to Cys in 40% of frames.

Since the single phosphate groups of pT-E-Y and T-E-pY tripeptides are tending to occupy the main active site, the corresponding secondary binding site is mostly empty, with the average distance between a center of mass of arginines ARG213-214 forming the second active site and phosphorylated residues of the tripeptides being about 10 Å (Supplementary Figure S10). The other reason for these long distances is that in most of the frames these arginines are not close to each other and do not form a tight active site. Residues interacting with only one arginine are located some distance from the center of mass.