Research ArticleIMMUNOLOGY

Robotic QM/MM-driven maturation of antibody combining sites

See allHide authors and affiliations

Science Advances  19 Oct 2016:
Vol. 2, no. 10, e1501695
DOI: 10.1126/sciadv.1501695


In vitro selection of antibodies from large repertoires of immunoglobulin (Ig) combining sites using combinatorial libraries is a powerful tool, with great potential for generating in vivo scavengers for toxins. However, addition of a maturation function is necessary to enable these selected antibodies to more closely mimic the full mammalian immune response. We approached this goal using quantum mechanics/molecular mechanics (QM/MM) calculations to achieve maturation in silico. We preselected A17, an Ig template, from a naïve library for its ability to disarm a toxic pesticide related to organophosphorus nerve agents. Virtual screening of 167,538 robotically generated mutants identified an optimum single point mutation, which experimentally boosted wild-type Ig scavenger performance by 170-fold. We validated the QM/MM predictions via kinetic analysis and crystal structures of mutant apo-A17 and covalently modified Ig, thereby identifying the displacement of one water molecule by an arginine as delivering this catalysis.

  • Immunology
  • antibodies


Immunoglobulins (Igs) play a key role in modern biopharmaceutics because of their unique ability to deliver the effective transfer of structural information to the target. The synergistic use of combinatorial and rational design has developed broad functionality within Igs for both binding and catalysis (15). This work has the potential to implement novel technologies of Ig engineering and expression (68). Approaches to Ig library screening using different anchors have become more efficient, with success for the development of an Ig binder or catalyst, especially in targeting toxic compounds (9). Nonetheless, a general restriction for “in vitro” immunization is its inability to effectively maturate Ig repertoires—a central feature of mammalian immunization that provides antibodies with optimum target specificity. We propose to resolve this deficiency in silico by predicting the binding/catalytic efficiencies using molecular dynamics (MD) (10) and quantum mechanics/molecular mechanics (QM/MM) approaches (9, 11), facilitated by supercomputing systems (12). Here, we develop a robotic approach to generate an artificially tailored prototype of a dead-end bioscavenger for an organophosphorus poison. This task has a tremendous social significance: Data from the World Health Organization identify an annual global mortality of >250,000 humans as a result of pesticide self-poisoning, which is 30% of all suicides (13, 14). Terrorism and illegal military action also sustain the need for defense against chemical weapons (15).


We chose as our initial target a wild-type Ig-paraoxonase (WTIgP) selected from a combinatorial library using a reactive phosphonate agent (fig. S1, A and B) (16, 17). Subsequently, WTIgP hydrolysis of paraoxon, a toxic organophosphate pesticide, was proven to proceed via covalent trapping (Fig. 1A and fig. S1C). We sought to enhance the reaction mediated by this naïve WTIgP by accelerating the first half-reaction (Fig. 1A). We thus needed to identify atomic details of the ground-state structure and covalent bond formation in this reaction.

Fig. 1 Mechanism of paraoxon hydrolysis by IgP.

(A) Paraoxon hydrolysis by WTIgP proceeds in two steps: covalent binding followed by dephosphorylation. ReAb, reactive antibody. (B) (Top) QM-computed mechanism of TS development for step 1 of WTIgP reaction with paraoxon: the classic SN2(P) mechanism with tbp geometry initiated by early proton transfer followed in 5 to 20 fs by LTS O–P bond formation. (Bottom) In the mutant reaction, Arg35 displaces one water molecule to coordinate strongly with the O==P of paraoxon; early proton transfer is followed in 10 to 35 fs by subsequent TS O–P bonding (H-bonds in red).

We first identified a candidate molecular mechanism of the antibody reaction using QM/MM procedures, focusing on the orientation of paraoxon in the active site of WTIgP. By molecular docking, we performed a cluster analysis of 2000 binding modes to achieve convergence of the mode population (fig. S2). The six top clusters were ranked in terms of binding energy (fig. S3): Clusters with energy of −6.3 to −5.8 kcal mol−1 showed stable positioning of the phosphate group, whereas those with energy of ≥−5.7 kcal mol−1 (fig. S3F) showed promiscuous orientation of phosphate. Although none of the six had the nucleophilic oxygen of Tyr37 positioned for in-line displacement (18) of the p-nitrophenolate (PNP) leaving group, conformational manipulation identified one cluster that could position paraoxon in a proreactive orientation. This adjustment called for rotation of only one P–O and one C–O bond relative to the fixed nitrobenzene ring (fig. S4).

A metadynamics QM/MM simulation at a semiempirical theory level showed reaction capable of proceeding from this higher-energy ground-state complex (GSC) via a trigonal bipyramidal (tbp) late transition state (LTS) to deliver the experimentally identified product (Fig. 1B). Moreover, the simulations identified direct proton transfer to the P==O group as an early transition state (ETS) component, a key feature of this phosphotriester reaction (Fig. 1B and movie S1). We next computed a four-dimensional (4D) energy plot, which has three collective variables (CVs) in the metadynamics simulation (Fig. 2C).

Fig. 2 The QM maturation scheme.

(A) Selected residues of amino acids that can theoretically donate a H-bond to paraoxon (light purple spheres), residues that can facilitate proton transfer at 5 Å from the Tyr37 hydroxyl group (red spheres). Positions capable of either of these two functions (yellow spheres). (B) Pyramid of mutant counts during in silico maturation. (C) Metadynamics simulation of paraoxon covalent binding pathway with WTIgP (left) and QM-matured Ser35Arg mutant (right). The metadynamics simulation used three CVs: proton distance to Tyr37 (TyrO---H), proton distance to phosphoryl oxygen (P==O---H), and breaking bond distance between leaving group oxygen and paraoxon phosphorus (P---O–PNP). A 4D plot with these three variables versus energy change is shown by color gradient (blue, low potential energy; red, high potential energy). Values for boxes (colored numbers) correspond to conformation populations at stages GSC, ETS, and LTS of the reaction (Fig. 1B). Marked decrease of conformation population in stage 2 in mutant form (58102 versus 9991) reflects reduced probability for the formation of H-bond between Tyr37 and phosphoryl oxygen.

We formulated a QM/MM maturation procedure to enable virtual screening of our mutant library, which focuses on four features, namely, (i) an unchanged mechanism of reaction, (ii) reprogramming of paraoxon positioning to form the most productive GSC, (iii) an effective proton transfer from Tyr37 to phosphate, and (iv) a cohort size manageable for in silico maturation. On the basis of the geometry of the hydrogen bonds formed, we searched for all virtual mutants in the active site, which deliver maximized contact between the substrate and amino acid residues of the active site (figs. S5 and S6). This identified 23 amino acids, which, when mutated to Arg, could form a H-bond (≤3.2 Å) with the phosphate moiety of paraoxon (Fig. 2A, light purple, red, and yellow spheres). We chose arginine for this search because its side chain gives the greatest inclusivity of amino acids for random selection. This design included all possible mutations that are capable of first shell polar interaction with the reacting species and gave a cohort size of 8.39 × 1029 (that is, 2023) combinations (Fig. 2B, purple segment). To reduce this to a more manageable size, we limited the selection of mutant amino acids to 11 polar residues and restricted permutations to four at a time. This cohort (1.7 × 109; Fig. 2B, blue segment, and Eq. 1) was still computationally unmanageable. We therefore programmed a smaller cohort of seven H-bond donor amino acids (Arg, His, Lys, Ser, Thr, Trp, and Tyr) at 19 chosen positions (Fig. 2A, light purple plus yellow spheres) and specified Glu, Asp, and Ser mutations to provide general acid-base catalysis for Tyr37 [Fig. 2, A (red plus yellow spheres) and B], because this catalysis is featured in many phosphoryl transfer reactions (19, 20). In combination, these mutations gave a manageable cohort of 167,538 candidates for development (Fig. 2B, cyan segment, and Eq. 2), using an MM Monte Carlo algorithm organized in the Rosetta package (21).

Our robotic in silico maturation process converted these virtual mutants into structural models, which were analyzed using a PyRosetta framework (21) and gave 375 hits (Fig. 2B, yellow segment). The first QM/MM iteration revealed that (i) systems with an anionic amino acid mutation did not deliver covalent reaction (Fig. 2A, red and yellow spheres) and (ii) successful runs resulted from cationic amino acid mutations involving residues Arg, Lys, and His in positions 35, 92, and 98 of the light chain and in positions 99 and 107 of the heavy chain. These residues were subselected for further individual QM/MM analysis (Fig. 2B, orange segment, and fig. S7) and MD simulation, which are capable of calculating the diffusion coefficient of paraoxon in the active center, thereby identifying Ser35Arg as the best mutant candidate for experimentation (Fig. 2B, red segment, and fig. S7A).

Expression and purification of the mutant Igs identified by the robotic maturation enabled us to evaluate them in vitro. Strikingly, steady-state kinetic analysis showed that Ser35Arg exhibits a 170-fold increase in k2/KD compared to WTIgP (229 M−1 min−1 versus 1.35 M−1 min−1; table S1). We mutated Ser35 to alanine and found that it was inactive to paraoxon. In addition, we found that Ser35Lys has the same k2 as Ser35Arg with an elevated KD, whereas Ser35His is 20-fold less reactive (table S1); both Ser35Lys and Ser35His were identified in silico to be reactive mutants. Pre–steady-state analysis of WTIgP, Ser35Arg, and Ser35His with paraoxon identified an induced fit activity for each protein, in agreement with the mechanism previously proposed for arylphosphonate binding (fig. S1A and table S2) (17). As a control, we analyzed the QM-inactive mutant Ser35Glu and found that it was experimentally inactive with paraoxon. It thus appears that an anionic residue inhibits the phosphoryl transfer reaction (table S1). We also found no turnover for the Ser35Arg and Ser35Lys mutations, showing that Ser35 is essential for the dephosphorylation of the covalent Ig intermediate and indicating that additional mutations are required to restore this turnover function. Ser35His is thus the only mutant that sustains paraoxon hydrolysis at the same rate as does WTIgP.

We determined two high-resolution crystal structures: one of Ser35Arg apolipoprotein (apo) antibody [Protein Data Bank (PDB): 5adp; 2.1 Å resolution] and one of Ser35Arg covalently bound to paraoxon (PDB: 5ado; 1.6 Å resolution) (fig. S8). The Arg35 residue is well defined, particularly in the paraoxon structure. Arg35 in the product structure forms an extended H-bonding network involving paraoxon, Asn105, and water molecules (Fig. 3A), as seen for WTIgP with an arylphosphonate ligand (PDB: 2xzc; Fig. 3B). This interaction evidently stabilizes the phosphorylated protein and narrows the cavity entrance. The phosphorylated Ser35Arg structure fully validates the selection criteria used for the QM/MM maturation analysis and its predictions. Moreover, superposition of the apo antibody and covalently modified antibody structures shows that Arg35 is displaced in providing space for the paraoxon substrate, with Cε translocating some 5.2 Å (Fig. 3C), a clear demonstration of induced fit for this antibody. Antibody catalysis can deploy conformational reorganization, as previously documented for Ig hapten binding (22), contradicting an historical criticism of this limitation for the viability of catalytic Igs (23).

Fig. 3 Structural analysis of mutant IgP.

(A) Active site and H-bonding at the active center of Ser35Arg (PDB: 5ado). Green, catalytic residues; gray, passive locating residues; black, diethyl phosphate; black dashes, H-bond network. (B) Active site and H-bonding at the active center of WT bound to a tropinyl phenylphosphonate derivative (PDB: 2xzc). Green, Tyr37; gray, passive locating residues; black, phenylphosphonate; black dashes, H-bond network. (C) Apo-Ser35Arg aligned with modified Ser35Arg (pink and magenta, with water molecules in pink), showing induced fit for Arg35 (Cε displaced 5.2 Å) (14).

We found that Ser35Arg-phosphate has an additional H-bond between the phosphoryl oxygen and the arginine guanidinium NH1, which should also be featured in the TS (Fig. 3A). Because there is no other H-bond interaction from this mutant to the Tyr-phosphate product (except a water molecule, also present in WTIgP, that coordinates both arginine NH1 and the P==O group), the 170-fold rate enhancement appears primarily as a result of a single H-bond from Arg35 that polarizes the P==O group, combined with the displacement of a water molecule. This point mutation to introduce an arginine has achieved the desired improvement in orientation of the substrate and electrophilic activation. QM/MM analysis enables a more detailed analysis of this result (Figs. 1B and 2C, fig. S7, and movie S2). This result, developed with minimal conceptual input following the selection of the maturation vector, has achieved a significant enhancement of the reactivity of the targeted phosphoryl group through robotic identification of a single point mutation. The fact that this result can be rationalized post hoc by crystallographic analysis is highly encouraging for the future development of this methodology.


Few enzymes catalyze phosphoryl transfer reactions using Tyr as an accepting moiety in covalent catalysis. They include topoisomerases and recombinases, whose active centers superpose remarkably well on our evolved Igs, whereas their overall molecular architecture differs significantly (Fig. 4). Our computer-matured Arg35 mutant maps much more closely onto the structure of these key enzymes than does WTIgP, but it does not share tyrosyl phosphate cleavage activity with them. Thus, the structures of topoisomerases and recombinases may hold catalytic features that will enable us to further mutate our engineered antibody to achieve complete turnover of organophosphorus agents.

Fig. 4 Comparison of structures of ligand-bound covalent complexes of Ser35Arg with topoisomerase and Cre-recombinase.

(A) Alignment of Ser35Arg-phosphate complex (PDB: 5ado; cyan) with topoisomerase IB (PDB: 2h7f; red). (B) Alignment of Ser35Arg-phosphate complex (PDB: 5ado; cyan) with Cre-recombinase (PDB: 1crx; pink). Catalytic residues were labeled with the same color code.

The in silico prediction of catalytic efficiency using supercomputer facilities has enabled us to optimize mutations in the Ig combining site. One of the main challenges in Ig-driven catalysis resides in its present inability to mimic the dynamics of highly evolved enzymes (24, 25), even allowing that this may not be an essential feature for many phosphoryl transfer enzymes (19). Here, we demonstrate that the QM/MM maturation approach not only is a valuable tool that enables the design of covalent capture of an organophosphorus agent by the antibody but also has the potential for stepwise extension to generate biocatalysis for substrate turnover. This maturation approach is capable of future development to more general applications, displaying wider diversity of Ig folds. The next supercomputer challenge will be tasked by simultaneous optimization of evolutionary deviations at the active center of the biocatalyst associated with the second sphere and longer-range interactions developed in the Ig folds.


Searching the best positioning of paraoxon in the active center of WTIgP

The AutoDock Vina package used the iterated local search global optimizer based on the Metropolis criteria for binding mode selection. It was natural to expect a result deviation between runs because the starting conformation was random. To overcome this uncertainty, we studied the saturation curve of the population of the best binding mode during increment program runs. Results are presented in fig. S2. The correct orientation of paraoxon for effective tyrosine attack is realized in approximately 5% of the analyzed variants. In its simplest form, the docked position of paraoxon appears to call for the unexpected “adjacent” attack of the tyrosine oxygen at phosphorus. It would deliver a pentacoordinate intermediate that can only bring the incoming and departing oxygens “in line” as a result of a pseudorotation (26). However, the conformational reorganization of C–O and O–P bonds of the lowest energy docking mode for paraoxon (shown in figs. S3A and S4) led to a higher-energy, reactive conformer by mobile repositioning of the diethoxy-phosphoryl group relative to the bound p-nitrophenolate moiety.

QM/MM calculations

The starting conformation WTIgP complex with paraoxon was taken from the docking studies described previously. The MM approach was applied to the MM part of the system using the parameters from parm99 force field with corrections (27). The QM subsystem was described with semiempirical Hamiltonian PM6-D3H4 (28). Side-chain atoms of Ser35 and Tyr37 (positions in light and heavy chains are denoted by numbers and italics, respectively) and the two closest water molecules were included in the QM part of the system. Each simulation system was filled with water molecules presented by the TIP3P model, and the total charge was neutralized with Na+ or Cl ions. Water and ions were equilibrated around the protein-paraoxon complexes by carrying out a 100-ps MD simulation, with a restrained position of WTIgP and paraoxon. The prepared systems were subjected to QM/MM simulation with the in-house modified GROMACS/MOPAC2012 package (29, 30). The time step used was 0.2 fs. Temperature coupling with Nosé-Hover scheme allowed observation of the behavior of systems at human body temperature. The total length of simulation was set to 20 ps. A metadynamics approach was used to overcome the activation barrier (31). In the first simulation run, we used one CV: the TyrO–H distance. After the formation of the TS, we used two CVs: distance P–O••••H and P••••O–PNP. Several metadynamics parameters (hill height, width, and addition frequency) were tested: A hill with a width of 0.01 nm and a height of 0.5 kJ/mol was chosen, and the hill addition frequency was selected as one hill every 50th step. To prevent paraoxon diffusion from attacking Tyr, this was distance-restrained to 2.5 Å with upper-wall CV as defined in PLUMED (32). This set of CVs led to approximately 50% probability of the TS formation.

We found that for the ETS (Fig. 1B), the approach of Tyr37 OH to phosphorus led, at an O–P distance of around 2.2 Å, to direct proton transfer from tyrosine to the phosphoryl oxygen. This is immediately followed (in ca. 25 fs) by the interaction of the tyrosine oxyanion with phosphorus to generate the LTS (Fig. 1B). Here, the developing O–P bond is around 1.9 Å, whereas the breaking bond has only stretched to 1.7 Å. Thus, this TS is very late, consistent with a relatively high-energy reaction that has only a single catalytic contribution.

Virtual library construction and screening

We estimated the number of virtual mutants as 1.7 × 109. The size of the library was calculated using Eq. 1Embedded Image(1)where k = (4 – number of sites occupied by polar amino acids) [amino acids (aa) = 11] and n = (23 – number of amino acid positions) that theoretically can donate a H-bond from the side chain to paraoxon.

We obtained 167,538 structural models of virtual mutants using MM Monte Carlo organized in the Rosetta package (21). The size of the library was calculated using Eq. 2Embedded Image(2)where n = (7 – number of combinatorial positions) for Glu, Asp, or Ser residues; r = (3 – number of sites) occupied by Glu, Asp, or Ser [amino acids (aa) = 3]; p = (7 – number of amino acid species) that theoretically can donate a H-bond from the side chain to paraoxon; and m = (19 – number of amino acid positions) that theoretically can donate a H-bond from the side chain to paraoxon.

Conformations of 167,538 possible mutants were analyzed using the PyRosetta framework (21), with a simple geometrical search, where side-chain torsion random rotation generates 100 conformations per mutant residue until the distance of the deviation of the conformation closest to the ligand was >0.2 Å. The top 50 conformations were subjected to geometry optimization, whereas the rest of the protein and ligand were fixed. The resulting 50 conformations were analyzed for the proximity of H-bond donor group to acceptors from the ligand (fig. S2). The best conformation in terms of contact count and with energy of >100 atomic units (estimated according to the Lennard-Jones equation; fig. S6 and Eq. 3) was checked for the possibility of forming a H-bond between the candidate mutant amino acid and paraoxon in terms of the distance between the two interacting partners. The resulting conformation was accepted if the H-bond distance was ≥3.2 Å. To select mutants for further QM/MM simulation, we introduced the following score function: score = N(contacts)/E(repulsion) × 100.

We applied a threshold score of >1.25 and selected 375 mutants.Embedded Image(3)

The selected 375 mutants were subjected to QM/MM metadynamics simulation as done previously for WTIgP, and the metadynamics parameters were chosen to simulate the TS formation with 50% probability for WTIgP.

These calculations are poorly parallelized because of the limits of the software used, and the best solution was to use a big cluster and also one computational core per reaction in a parallel manner. From this point of calculation, all QM/MM runs were performed using the Lomonosov supercomputer with 52,168 cores, which has a 30 times overhead in power needed.

After a major run, the top 10 candidates were selected for precision evaluation. Here, we performed 100 runs for each mutant to achieve convergent results and better sampling. In this case, upper-wall CV for the distance between phosphorus and Tyr-oxygen was increased to 2.75 Å. This value corresponds to 20% probability of the TS formation for WTIgP.

Energy landscape calculation

3D representation of the energy landscape of reaction for mutants with paraoxon as a function of three CVs was built. Because our reaction actually involves two steps, we average here two independent runs with different sets of CVs. We reconstructed the energy plot on the basis of a histogram occurrence of states, F(s) = −kBT log[P(s)], where P is the probability to meet an exact set of CVs. Each conformation from a metadynamics run has a set of values of CVs. We selected three intervals of CV values corresponding to reaction stages from Michaelis complexes and from activated Michaelis complexes to the TS. We naturally assigned and counted conformations on the basis of CV values for the different stages. For example, the basic Michaelis complex corresponds to CV(Tyr–O…H) < 1.2 Å, CV(P==O…H) > 2.5 Å, and CV(P…O–PNP) < 1.9 Å. The rotationally activated Michaelis complex corresponds to CV(Tyr–O…H) < 1.2 Å, CV(P==O…H) < 2.5 Å, and CV(P…O–PNP) < 1.9 Å. The transition complex corresponds to CV(Tyr–O…H) > 1.2 Å, CV(P==O…H) < 1.2 Å, and CV(P…O–PNP) < 1.9 Å.

Movie making

Video files were generated on the basis of QM metadynamics trajectory using the PyMOL software.

Calculation of diffusion coefficient

Structure modeling. Paraoxon coordinates were built with Open Babel from SMILES notation. Molecule geometry was optimized at a B3LYP 6-31G(3d,2f) level followed by point atomic charge calculation. Point charges were derived from restrained electrostatic potential (RESP) calculated on the same level of theory with R.E.D. (RESP and ESP charge Derive) utility (33). Molecular topology with optimized coordinates and derived atomic point charges was built by ACPYPE (34). Initial structures of all complexes were built on the basis of the position of paraoxon in the rotationally activated state of the Michaelis complex. All simulation and analysis studies were performed using the GROMACS 5.0 software package.

Classical MD simulations. Explicit solvent simulations in parm99sb-ILDN force field (35) were performed at T = 300 K under control of a velocity-rescaling thermostat (36), with isotropic constant-pressure boundary conditions under the control of the Berendsen algorithm of pressure coupling (37), and with application of the particle mesh Ewald method (36) for long-range electrostatic interactions. A triclinic box of the TIP4P (36) water molecules was added around the protein to a depth of 15 Å on each side of the solute. Charges of systems were neutralized by addition of chloride anions, resulting in ~0.15 M chloride ions. We used two temperature coupling groups in the simulations. The first was a modified protein-paraoxon complex, and the second was water with anions. Analyses were carried out using the tools from the GROMACS package with Python scripts. Diffusion coefficients were calculated on the basis of 1-ns-long trajectory with the Einstein relation.

Construction of the expression vectors and site-directed mutagenesis

Recombinant Fab fragment of WTIgP and its mutants were produced in the methylotrophic yeast Pichia pastoris GS115 (Invitrogen) using the modified expression vector pPICZα/Jk1 (38) based on pPICZα (Invitrogen). Mutations at light chain Ser35Arg, Ser35His, Ser35Glu, Ser35Lys, Ser35Ala, and Pro98Arg were created by an overlap extension with two DNA fragments obtained by polymerase chain reaction (PCR), using the pPicZα/A17/Lch plasmid, which contained A17 (WTIgP) κ light chain, as a template and the following primers: a/c and b/d for the pair encoding Ser35Arg mutation, a/f and b/e for the pair encoding Ser35His mutation, a/h and b/g for the pair encoding Ser35Glu mutation, a/i and b/j for the pair encoding Ser35Ala mutation, a/k and b/l for the pair encoding Ser35Lys mutation, and a/m and b/n for the pair encoding Pro98Arg mutation. Overlap extension amplification for each mutation was completed using the primer pair a/b. Final PCR products for light chain mutants were digested by BspMI and XhoI restriction endonucleases and used to replace the insert between BsmBI and XhoI sites in the pPICZα/Jk1/κ. Mutation in heavy chain Leu99Arg was created by an overlap extension with two DNA fragments obtained by PCR, using the pPicZα/A17/Hch plasmid (WTIgP heavy chain containing the CH1 constant domain with hinge region) as a template and the primer pairs o/p and q/r. Overlap extension amplification was completed using the primer pair o/q. Final PCR product for heavy chain mutant was used to replace the insert between KpnI and PstI sites in the pPicZα/A17/Hch plasmid. For rapid selection of mutated clones, the restriction site of KpnI was deleted in the sequence in the case of Ser35Arg mutation. The restriction sites of Eco105I and ApalI were inserted in the sequence in the case of Ser35Glu and Ser35His mutations, respectively. The restriction site of SspI was inserted in the sequence in the case of Ser35Ala. The restriction site of SnaBI was inserted in the sequence in the case of Ser35Lys. The restriction site of Pvu II was inserted in the sequence in the case of Leu99Arg mutation. Last, the restriction site of Acc65I was inserted in the sequence in the case of Pro98Arg mutation. All constructs were verified by DNA sequencing.

The additional primer sequences were as follows: (a) ttgatctcgagcttggtccctcc, (b) gaaggacctgcaggaaagacagtctgtgctgacgcag, (c) ttatgtacgctggtatcagcagc, (d) gctgctgataccagcgtacataa, (e) tgctggtaccagtgcacataattattc, (f) aataattatgtgcactggtaccagcag, (g) gctgctggtaccattctacgtaattattc, (h) aataattacgtagaatggtaccagcagc, (i) aataattatgtagcttggtaccagcagctc, (j) accaagctacataattattcccaatattggagct, (k) aataattacgtaaagtggtatcagcag, (l) ctgataccactttacgtaattattccc, (m) aatcgcgtattcggcggaggtaccaag, (n) cttggacctccgccgaatacgcgatt, (o) aaaagacaggtgcagctgcaggag, (p) cgtgcgtcctgcacagtaatacacagctgt, (q) gcctgagttccacgacaccgtca, and (r) actgtgcaggacgcacgcagtc.

Protein expression and purification

Procedures of electrocompetent cell preparation, electroporation of P. pastoris GS115 cells, Mut+ or Muts phenotype determination, and selection on Zeocin followed Invitrogen protocols. Analytical or large-scale expression of recombinant WT and its mutants was performed in cultures of BMGY and BMMY media according to Invitrogen protocols. Methanol was added every 24 hours after induction (up to 0.5%). WTIgP and its mutants were purified as described previously (38). The purity and identity of the eluted Fabs were tested by 12% SDS–polyacrylamide gel electrophoresis with Coomassie staining and Western blot analysis. Horseradish peroxidase–conjugated anti-FLAG and anti-human κ light chain antibodies (Sigma) were used to detect heavy chain constant fragment (CH) and light chain constant fragment (CL), respectively.

Evaluation of kinetic and thermodynamic parameters for the WTIgP and its mutant reactions

Kinetic measurements were made as described previously (16, 17). Briefly, reactions of WT and mutants (3 to 10 μM) with phosphonate X or paraoxon over a concentration range of 10 to 500 μM were carried out in 0.1 M sodium phosphate buffer (pH 7.4) at various temperatures. Reaction rates were determined from the changes in absorbance at 400 nm that resulted from PNP formation, and the rate constants were calculated using a PNP extinction coefficient ε of 12,300 M−1 cm−1. Modifications of rate constants k1 were estimated by Kitz-Wilson analysis (39) using the DynaFit software (40).

Stopped-flow measurements with fluorescence detection were carried out using an SX.18MV stopped-flow spectrometer (Applied Photophysics). The dead time of the instrument is 1.1 ms. The antibody concentration was fixed at 10 μM, and that of paraoxon was varied from 5 to 1000 μM. The reaction buffer contained 50 mM sodium phosphate (pH 7.4), 150 mM NaCl, and 40% glycerol at 1oC. Glycerol was used to retard the reaction rate according to Pocker and Janjic (41). The fluorescence of the reaction mixture was detected using λex = 290 nm (bandwidth, 5 nm) and emission filter with a cutoff of 320 nm. Each trace shown is the average of four individual experiments.

Accurate modeling of all kinetic traces is obtained by simulation with the DynaFit software, using an appropriate reaction mechanism with a single set of constants. The conformity of the fitted kinetic models to the experimental data is assessed by monitoring residuals against time for different scheme fits.

The integral activation energy for the phosphotyrosine intermediate formation was calculated as 21 kcal/mol. This estimate looked very realistic and was comparable with the estimate for the noncatalytic process. All predicted mutants of WTIgP were engineered, expressed, and subjected to a detailed study by steady-state (table S1) and pre–steady-state analyses (table S2). The very close “experimental” and QM/MM-estimated energetic profiles for WTIgP provide a good quantitative outturn for our developed approach. Because of simplifications necessarily implemented in the simulation methods, it was satisfactory to find the best variant by experimental screening from the narrow repertoire of mutants.

Mass spectrometry

To prove the covalent phosphotyrosine intermediate formation for Ser35Arg, we performed mass spectrometric analysis of the tryptic hydrolysate of the Ser35Arg mutant before and after interaction with paraoxon. The Ser35Arg antibody (10 μM) was incubated with paraoxon [2 mM in 0.1 M phosphate buffer (pH 7.4)] at 37°C for 16 hours. Samples were treated with sequencing-grade trypsin according to the Bruker applications guide and analyzed directly using UltraFlex II MALDI-TOF/TOF mass spectrometer (Bruker Daltonics). Analysis was performed with flexAnalysis 3.3 (Bruker Daltonics). The mass difference (M/e, 136) of the multicharged peaks corresponded to the mass of a diethoxyphosphoryl residue (fig. S9).

Crystallization and structural determination of Ser35Arg Fabs, unmodified and modified by paraoxon variants

Crystals of both the Ser35Arg variants, unmodified and modified by paraoxon, were grown using the hanging-drop vapor diffusion method by mixing equal volumes of the protein [8 mg/ml in 50 mM tris-HCl buffer (pH 7.4)] and precipitant solution [0.15 to 0.2 M MgCl2 and 12 to 18% (w/v) polyethylene glycol (PEG) 6000 in 0.1 M N-(2-acetamido)iminodiacetic acid (pH 7.5)]. The crystallization conditions for these variants were very similar to the ones for the native WTIgP protein (17). Rod-shaped crystals of approximate dimensions (0.4 mm × 0.2 mm × 0.2 mm) were obtained after few days. For the unmodified protein, x-ray data were collected on the European Molecular Biology Laboratory (EMBL) beamline X13 at the DORIS storage using a 165-mm MARCCD detector at a wavelength of 0.8123 Å. Images of 0.45° oscillation were collected for a total of 150° of rotation, using 45-s exposure per image. Data with the paraoxon-modified variant were collected at P14 at the PETRA III storage ring [German Synchrotron Research Center (DESY)] at a wavelength of 0.8266 Å, using a Pilatus 6M detector. The chosen oscillation range was 0.1° for a total of 360° of rotation, using 0.2-s exposure per image. Data were collected at cryogenic temperature (100 K), and the mother liquor solution, supplemented with 20% (v/v) PEG 400, was used as a cryoprotectant. Diffraction data were indexed and integrated with XDS (42) and scaled using SCALA (43). Values of I/σ(I) and the CC1/2 were used as a guide to determine the resolution cutoff for the paraoxon-modified Ser35Arg antibody, whereas we used data of the highest collected resolution for the unmodified antibody (tables S3 and S4) (44). The coordinates and experimental structure factors of the x-ray structures are deposited in the Research Collaboratory for Structural Bioinformatics PDB, with code 5adp for the unmodified Ser35Arg and code 5ado for the paraoxon-modified variant.

Both structures were isomorphous with the WTIgP crystal structures, belonging to the P1 21 1 space group. The asymmetric unit of the crystal contained one Fab molecule and four domains (VL, CL, VH, and CH1), which have the canonical β-sandwich Ig fold (fig. S5). The Ser35Arg structures were solved by molecular replacement using MOLREP (45), with the WTIgP structure [PDB: 2xza (17)] as a search model. The refinement was carried out with PHENIX (46), including TLS (Translation-Libration-Screw-rotation model) refinement. Solvent molecules were located automatically using PHENIX, confirmed by visual inspection, and refined with unit occupancy. For the unmodified Ser35Arg, the final model comprises residues 2 to 223 of the heavy chain and residues 2 to 218 of the light chain, with the exception of residues 101 to 106 of the H-CDR3 (heavy chain complementarity-determining region 3) loop and residues 135 to 139 of the heavy chain. For the paraoxon-modified variant, the final model consisted of residues 2 to 222 of the heavy chain and residues 4 to 216 of the light chain, omitting residues 101 to 104 of the H-CDR3 loop and residues 134 to 139 of the heavy chain. A stereochemical analysis of the structure using PROCHECK (47) showed that 98.1 and 98.9% of the residues were in the most favored regions of the Ramachandran plot, and 1.9 and 1.1% were in favored regions for the unmodified and modified structures, respectively. The refinement was converged at Rfactor of 17.9% and Rfree of 22.0% for the unmodified structure and at Rfactor of 16.8% and Rfree of 19.9% for the paraoxon-modified variant.

As expected, the overall structures of the Ser35Arg mutant share high similarity with the WTIgP [root mean square deviation (RMSD) of 0.4 Å when aligned using the LSQKAB algorithm)]. The light chain is practically identical in all the structures (RMSD of about 0.2 and 0.3 Å for Ser35Arg unmodified and Ser35Arg modified with paraoxon, respectively). There is a higher deviation for the heavy chain (RMSD of about 0.5 Å when aligned using the LSQKAB algorithm), and this difference is mainly due to the fact that, in the Ser35Arg structures, the CDR loops of the heavy chain, and especially the H-CDR3 and H-CDR1, are highly disordered in contrast to the WTIgP structures (table S5). In contrast to this, the CDR loops of the light chain of Ser35Arg have similar normalized atomic displacement parameter values with the WTIgP structures.

Structural comparison of Ser35Arg mutant with topoisomerases and recombinases

Alignment of the x-ray structures of ligand covalent complexes of the paraoxon-modified Ser35Arg antibody with topoisomerase IB (PDB: 2h7f) and Cre-recombinase (PDB: 1crx) was performed using the “pair fitting” function of the PyMOL software. Fitting was performed using five atoms: Tyr37-oxygen, P (phosphorus atom for paraoxon or DNA phosphate group), and O1, O2, and O3 (phosphoryl oxygens for paraoxon or the DNA phosphate group).


Supplementary material for this article is available at

fig. S1. Chemical structures of compounds used in the study.

fig. S2. Saturation curve of population of the best binding mode during incremental program runs.

fig. S3. Cluster analysis of paraoxon binding modes in the WTIgP active center [PDB: 2xza (17)].

fig. S4. Development of Michaelis conformation of paraoxon for in-line attack by Tyr37.

fig. S5. Typical set of 50 optimized conformations for the Ser35Arg mutant.

fig. S6. Graphical estimation of threshold for optimal conformation selection.

fig. S7. QM/MM analysis of top 9 mutants.

fig. S8. Overall structure of the Ser35Arg unmodified antibody (A) (PDB: 5adp) and paraoxon-modified antibody (B) (PDB: 5ado) structures.

fig. S9. Mass spectrometric analysis of the tryptic hydrolysate of Ser35Arg (WYQQLPGTAPK) antibody interaction with paraoxon.

table S1. Kinetic parameters for the interaction of paraoxon and phosphate X with WTIgP and seven mutants.

table S2. Pre–steady-state analysis of WT, Ser35Arg, and Ser35His reaction with paraoxon.

table S3. Data collection and refinement statistics for the unmodified Ser35Arg apo antibody.

table S4. Data collection and refinement statistics for paraoxon-modified Ser35Arg antibody.

table S5. Comparison of atomic displacement parameter and normalized atomic displacement parameter values for the CDR regions of the WTIgP and Ser35Arg reactibody structures.

movie S1. QM-computed mechanism of TS development for step 1 of WTIgP reaction with paraoxon: the classic SN2(P) mechanism with tbp geometry initiated by early proton transfer followed in 5 to 20 fs by LTS O–P bond formation.

movie S2. QM-computed mechanism of TS development for step 1 of the mutant Ser35Arg reaction with paraoxon: Arg35 displaces one water molecule to coordinate strongly with O==P of paraoxon; early proton transfer followed in 10 to 35 fs by subsequent TS O–P bonding (H-bonds in red).

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We thank G. Bourenkov and J. Kallio for technical assistance and A. Zalevsky for GROMACS/MOPAC2012 interface development. Funding: This work was supported in part by a German-Russian collaborative grant [ID RFMEFI61614X0009/BMBF (Federal Ministry of Education and Research) project: 05K14YEB] for the groups of M.W. and A.G.G. X-ray data were collected on the EMBL/DESY beamlines P14 (Petra III) and X13 (DORIS III), Hamburg, Germany. Author contributions: G.M.B., A.G.G., and R.A.L. conceived the study. G.M.B., S.D.C., A.G.G., A.V.G., R.A.L., I.V.S., and M.W. planned the experiments, interpreted the results, and wrote the manuscript. M.W. and S.D.C. planned and performed the x-ray crystallographic analysis. A.V.G. designed and executed the computational work described. A.A.B., S.D.C., A.V.G., I.N.K., N.A.P., I.V.S., A.V.S., Y.P., and O.I.Z. carried out the experiments. Competing interests: The authors declare that they have no competing interests. Institutional review board statement: All experimentation with hazardous chemicals was in accordance with the institutional guidelines. Data and materials availability: Data for the protein structures have been lodged in the PDB under codes 5adp for the mutant A17 S35R apo antibody and 5ado for the WT antibody modified by paraoxon. Additional computational information is available from A.V.G. upon request.

Stay Connected to Science Advances

Navigate This Article