Deciphering the Structural Effects of Activating EGFR Somatic Mutations with Molecular Dynamics Simulation

Molecular dynamic simulation is a computational technique that can be used to explore molecular motions and can reveal conformational changes crucial to understanding biochemical and cellular function. Probing the range of dynamic motions accessible to a macromolecule is difficult. Using the results from molecular dynamic simulations in conjunction with experimental data allows assessment of their functional significance.

We will demonstrate using molecular dynamic simulations how mutations observed in cancer patients affect EGFR targeting confirmations and ligand binding. To prepare the wild-type apo active EGFR our kinase structure, open the Kymera visualization program, and under the file menu, click fetch by ID.Select the Protein Data Bank database and specify the Protein Data Bank 2GS2 code. To build the missing structural 2GS2 elements, acquire these segments from other EGFR structures.

To construct the five residue glutamate 746 to alanine 750 ELREA deletion in EGFR, click favorite, sequence, and show sequence to open the wild-type 2GS2 sequence. And in the resulting sequence window, click edit and add sequence to select the deletion mutant FASTA format sequence. In the alignment window, select the structure and model or homology.

In the popup window, specify the 2GS2 composite structure as the template and the mutant sequence as the query to be modeled, then select a mutant model from the resulting models based on the zDOPE score and visual inspection. To prepare the wild-type apo inactive EGFR kinase structure, open Protein Data Bank structure 2GS7 and add the missing segments from the other EGFR structures, modeling the deletion mutant form as demonstrated. To prepare the ATP-bound wild-type active EGFR kinase structure, use Protein Data Bank structure 2ITX as the principle structure, building the missing segments using other EGFR structures and modeling the deletion mutant form using the modeler as demonstrated.

To construct the wild-type EGFR asymmetric dimer structure, open 2GS2 in Kymera and click tools, higher order structure, and unit cell to convert the structure to the biological assembly that contains the activator and receiver kinases in the asymmetric arrangement. Select the 2GS2 structure and enter make copies, then select and save a single asymmetric dimer from the multiple copies of the dimer resulting from the symmetry operations. To build the alanine 702 valine mutant, select tools, structure editing, and rotamers to replace alanine 702 with valine.

Open the structures in Maestro and click the protein preparation wizard button, then select add hydrogen atoms and fill in missing side chain atoms and click pre-process. To determine the protonation states of ionizable residues at pH 7.0, click refine and use PROPKA to optimize the orientation of asparagine, glutamine, and histamine residues for hydrogen bonding, then minimize the structure. To set up the simulation system, open the LEAP program and import the Amber FF14SB force field and TIP3P water molecules.

For the ATP-bound systems, import parameters for ATP and load the structure. Solvate the structure in an octahedral box with explicit TIP3P water molecules that extends 10 angstroms in all directions from the surface atoms of the protein. Check the belt system and add the necessary ions to neutralize it.

To sufficiently model biomolecular systems, add additional sodium and chloride atoms to the simulation box to bring the system salt concentration to 0.15 molar, then generate and save the topology and coordinate files of the system to serve as inputs for the subsequent production simulation. Using Amber, initially energy minimize the simulation system to circumvent any unfavorable configurations. In the minimization input file, adjust the maximum cycle variable for the total minimization cycle and the number of cycles to indicate the number of cycles for the steepest descent algorithm.

Use the restraint weight variable to apply the restraint force on the solute atoms specified by the restraint mask parameter. Carry out the minimization in multiple steps, gradually lowering the restraint applied on the solute atoms from 25 to 0 kilocalories per mole angstrom square, then use the command to run the minimization. Heat the system for 100 picoseconds from 0 to 300 Kelvin and use the commands to set a 10 kilocalorie molar per square angstrom restraint on solute atoms.

Then use the command to carry out the heating. Equilibrate the system for 900 picoseconds under an isothermal isobaric ensemble and set a nine angstrom distance cutoff for long range electrostatic interactions, gradually lower the solute atom restraint to 0.1 kilocalorie per mole angstrom square and finalize the equilibration with an unrestrained five nanosecond simulation. Run the equalization with the command as indicated.

Adjust the text to allow the production simulation to be carried out for 100 nanoseconds with the conformations being saved every 10 picoseconds, then run the simulation with the command as indicated. To visualize the conformation sample during the wild-type and mutant EGFR kinase simulations, open the Amber topology files and the corresponding trajectory files in Visual Molecular Dynamics. And using convenient secondary structure representations, analyze the overall structural dynamics of the proteins from the recorded trajectory.

Then view specific interactions between atoms and residues of interest, such as the catalytically essential lysine 745-glutamate 762 salt bridge. Alternatively, save multiple conformations sampled during the simulation in Protein Data Bank format and open the conformations in Kymera. Use the matchmaker option to superimpose the structures on the initial or median structure and display the median structure in solid and the rest of the aligned structures in faded white to allow visualization of the recorded structural movements with more clarity.

To analyze the global stability of the wild-type and mutant EGFRs and to examine the flexibility of the different structural units, import the Amber typology and corresponding trajectory files. In the root mean square deviation input file, indicate the backbone atoms of the initial structure as reference for the root mean square fitting. In the root mean square fluctuation input file, indicate the C-alpha atoms of the initial structure as reference for the root mean square fitting.

Then run the analysis with the CPPTRAJ program and plot the output data. Alternatively, to align the conformational ensembles and to color each residue based on the C-alpha atom root mean square deviation, open the conformations in Kymera and align them with the matchmaker option. Select tools, depiction, and render by attribute.

Select residues of the conformational ensembles and set the C-alpha root mean square deviation as the attributes, then click OK.The chain trace of the conformations will be colored blue, white, or red reflecting the regions of high, medium, and low structural stability respectively. To analyze the hydrogen bond interactions between ATP and both wild-type and deletion EGFRs, prepare a CPPTRAJ script to carry out this task. Specify the analysis for the intermolecular hydrogen bonds only with the nointramol variable and define a hydrogen bond with a donor acceptor distance of less than or equal to 3.5 angstroms and a bond angle of greater than or equal to 135 degrees.

To assess the intramolecular interactions, for example, between the catalytically important lysine and glutamate residues, specify lysine as the hydrogen donor and glutamate as the acceptor residues and run the script as indicated to allow analysis of the result. To compute the binding free energies between ATP and both wild-type and deletion EGFRs, prepare the gas phase ligand receptor and ligand receptor complex Protein Data Bank files in the LEAP program and set ATP as the ligand and EGFR as the receptor. Set the generalized born radio value to M Bond I2.Then generate the Amber topology and coordinate files for the gas phase Protein Data Bank files.

Similarly, to calculate binding free energies between the activator and receiver kinases of wild-type and alanine 702 valine EGFRs, specify the receiver kinase as the ligand and the activator kinase as the receptor and save the corresponding topology and coordinate files. Prepare a molecular mechanics generalized born surface area input file and set the IGB value to two and the saltcon to 0.1. Then using the mmpbsa.

py script available in Amber, enter the command as indicated to execute the binding energy calculations and analyze the output data. During this representative 100 nanosecond simulation, the alanine 702 valine mutant showed increased conformational stability of the juxtamembrane B segment likely due to tighter hydrophobic interactions as compared to the wild-type EGFR. The alanine 702 valine mutant also exhibited a lower free energy of binding between the activator and receiver kinases relative to the wild-type EGFR, representing more favorable dimer interactions that maintain the active EGFR kinase conformation.

Simulation of the deletion mutation resulted in lower C-alpha atom fluctuations of the functionally key alpha C-helix as compared to the wild-type EGFR, prolonging the time of the EGFR active state. The deletion mutation also resulted in frequent formation of hydrogen bonds between the side chain polar atoms of lysine 745 and glutamate 762, a key interaction for EGFR enzymatic activity compared to wild-type EGFR. Additionally, the number of hydrogen bonds between ATP and EGFR were greater for the deletion mutant than for the wild-type EGFR.

The deletion mutation also resulted in an inward movement of the EGFR inactive state alpha C-helix, a structural change expected during the transition to the active state. In contrast, the alpha C-helix of wild-type inactive EGFR maintained its initial conformation. To assess the impact of the mutations, it is critical to select the appropriate conformational states of the structures and to properly prepare and equilibrate these structures.

It is important to combine molecular dynamic simulations with experimental studies as the synergy between these techniques benefit the interpretation of the results and can inform additional wet lab experiments.