Computational protocol: De novo active sites for resurrected Precambrian enzymes

Similar protocols

Protocol publication

[…] To complement the experimental data, we also performed MD simulations on the wild-type B. licheniformis, E. cloacea (NMC-A), P. vulgaris (Bla-B), TEM-1, ENCA, GNCAMP, GNCA4 and PNCA β-lactamases, as well as the corresponding GNCAMP-W229D, GNCA4-W229D/F290W and PNCA-W229D variants (PDB IDs: 4BLM, 1BTL, 3ZDJ, 4B88, 5FQQ, 4C6Y, 4UHU, 5FQI respectively). We also performed simulations on the BL-W2299D, NMC-A-W229D, BLA-B-W229D, TEM-1-M182T/W229D and PNCA-W229D variants; however, as crystal structures of three enzymes were unavailable, we manually inserted these substitutions into the wild-type enzyme using PyMOL’s Mutagenesis Wizard. The specific GNCA variants were selected on the basis of their observed Kemp eliminase activities (), and the variant simulations were performed both in complex with the transition-state analogue (TSA) 5(6)-nitrobenzotriazole (PDB IDs for the β-lactamase-TSA complexes: 5FQJ for GNCAMP-W229D, 5FQK for GNCA4-W229D/F290W), as well as in the TSA free form. For all other variants, the TSA was manually placed in the newly created cavity in agreement with structural alignment of the crystal structures from the TSA bound ancestral variants. All MD simulations were performed with the GPU implementation of the AMBER16 simulation package, and the ff14SB force field. The Sander module was used for the initial minimization and equilibration runs, and PMEMD for the production runs in explicit solvent. To obtain ff14SB compatible parameters for the TSA, we parameterized this compound using ANTECHAMBER and Gaussian09 (ref. ). Partial charges were obtained using the standard RESP fitting procedure, and the ligand atoms were described by the GAFF force field. All non-standard parameters generated by the GAFF force field are provided in the .All model systems were generated with LEaP. Each model system was placed in an octahedral TIP3P water box, which allowed the box to extend to at least 11 Å from the solute in each direction. The protonation state of all ionizable residues was determined by empirical pKa calculations using PROPKA 3.1 (ref. ), and by visual inspection (in particular in the case of histidine side chains). On the basis of this, ionizable residues were all kept in their standard protonation state at physiological pH, and all histidine residues were kept neutral. This led to total system charges of −7, −5, −7 and −8 for TEM-1, BL, GNCA4 and GNCAMP respectively. These systems were then neutralized by addition of the appropriate number of Na+ counterions depending on the total charge of the system. The solvated systems were then subjected to a two-step minimization procedure to remove clashes between the water molecules and the solute. This was comprising of 50 steps of steepest descent and 200 steps of conjugate gradient minimization using 25 kcal mol−1 Å−2 harmonic positional restraints on all solute atoms, followed by 50 steps of steepest descent and 200 steps of conjugate gradient minimization with weaker 5 kcal mol−1 Å−2 harmonic positional restraints on all solute atoms. After this initial minimization, the following four sequential equilibration steps were performed: a 50 ps NVT simulation to increase the thermostat target temperature using the Berendsen thermostat and pressure control algorithms with 0.5 ps time constants for both the bath coupling and pressure relaxation; a 50 ps NPT simulation at a constant isotropic pressure of 1 atm to adjust the density of the system to 1 g cm−3; five 50 ps NVT simulations in which remaining the 5 kcal mol−1 Å−2 harmonic positional restraints were progressively decreased in 1 kcal mol−1 Å−2 increments, allowing us to finally perform a 50 ps NVT simulation without any restraints on the system.For each system, we performed three independent equilibrations using different initial random seeds, the starting points for which were obtained by increasing the target thermostat temperature from 100 to 299.9 K, 300.0 and 300.1 K, respectively, to obtain three independent simulations. This was followed by performing 100 ns production simulations for each of the three replicas. All simulations were performed using a 2 fs time step, saving snapshots every 20 ps. The SHAKE algorithm was used during the dynamics to constrain all bonds involving hydrogen atoms. Short-range non-bonded interactions were calculated subject to an 8 Å cut-off radius. Long-range interactions were described using the particle mesh Ewald method. The temperature was kept constant at 300 K using Berendsen’s weak coupling algorithm. All production simulations were performed using NVT conditions. Subsequent analyses of the MD trajectories were performed using CPPTRAJ. The first 40 ns of each simulation was discarded from the analysis as equilibration time, and all RMSF and clustering data shown in the paper are obtained as averages over the last 60 ns of three independent equilibration runs for each system (that is, 180 ns total simulation time per system). The backbone r.m.s.d. plots shown in demonstrate that all systems remained stable over the remaining simulation time. In , we compare the average RMSF values per residue with the corresponding averages of the B-factors derived from X-ray crystallography. Although B-factors have been often interpreted in dynamical terms, they are affected by a large number of other factors, including the resolution limit, radiation damage, crystal lattice defects, rigid body motions, occupancy levels and refinement artifacts. For the protein systems studied here, we observe an approximate congruence between the B-factors and RMSF values that appears comparable to the congruence reported for other protein systems in the literature. […]

Pipeline specifications

Software tools PyMOL, AMBER, PROPKA
Applications Small-angle scattering, Protein structure analysis
Chemicals Amino Acids