Atomistic simulations reveal impacts of missense mutations on the structure and function of SynGAP1

In silico mutagenesis and molecular modelling

The AlphaFold model for SynGAP α2 isoform-1 (UniProt: Q96PV0–01; AF-Q96PV0-F1-model_v3.pdb; https://alphafold.ebi.ac.uk/; 22 January 2022) [28]) was edited to contain reliable (pLDDT ≥50) N-terminal region and the Gly-rich loop (res. 198–730; Fig. 1B and C). Altogether 204 missense variants from ClinVar (https://www.ncbi.nlm.nih.gov/clinvar/; 21 April 2023) [30] and seven from other sources (e.g. MASTERMIND; https://mastermind.genomenon.com/; 21 April 2023) [31] were generated using MODELLER10.3 [32] (Table S1). SynGAP-membrane positioning (Fig. 2A and Fig. S2A) was acquired by aligning the Type II C2 domain against a Protein Data Bank-entry (PDB; https://www.rcsb.org/) [33] with settled C2 domain-membrane orientation at the Orientations of Proteins in Membranes (OPM; https://opm.phar.umich.edu/; 11/10/2022; PDB: 4RJ9) [34, 35] using VERTAA in BODIL [36] (Fig. S2B–C available online at http://bib.oxfordjournals.org/). SynGAP-RasGTPase complex was modelled by aligning the GAP domain with a human GAP-RasGTPase structure (PDB: 1WQ1; Fig. 3A and B) [37], and, at the active site, the GDP-AlF₃ was converted to GTP in MAESTRO (Schrödinger Release 2022–3: Maestro, Schrödinger, LLC, New York, NY, 2022). A more detailed description is given in the Supplementary Information (SI).

Molecular dynamics simulations and protein stability calculations

SynGAP-solvent systems were prepared for the WT protein (Fig. 1B) and missense variants with TLEAP in AMBERTOOLS22 [38] and simulated using AMBER22 [39], each running for 3 x 150 ns. The AMBERff19SB and OPC force fields were used for the protein and solvent, respectively [40, 41]. The SynGAP-membrane (Fig. S2A), SynGAP-RasGTPase-membrane complex (Fig. 3A) and intracellular membrane only [42, 43] systems were modelled using CHARMM-GUI and MD simulated 130, 300, and 150 ns, respectively, using GROMACS2021 [44] and CHARMM36ff [45]. Ras was Cys-palmitoylated at residues 181 and 184 (Figs. 3A and Fig. S3). Protein folding stability calculations for the variants were performed using FoldX5.0 [46] based on frames from the 150-ns WT-SynGAP-solvent simulations [47] and, likewise, using Rosetta Cartesian ΔΔG protocol [48–54]. Protocols (Fig. S4) and input files (supplementary_files.zip) are given in the SI.

Analysis and figure preparation

The MD trajectories were examined, and figures were generated using MOL^*4.0.1 [55], VMD1.9.4a12 [56], PYMOL2.3 [57], and MAESTRO. The trajectory analysis, including Root-Mean-Square-Deviation (RMSD; Figs. S5–S6) and atomic fluctuation (Fig. S7) calculations [58], was performed using CPPTRAJ in AMBERTOOLS22 [38]. With SynGAP-membrane simulations, the lipid bilayer metrics were analysed using GROMACS2021 [44] tools (Fig. S8).

Resilience of the SynGAP model in simulations

The SynGAP-solvent MD simulations, probing the structural effects of the missense mutations (N = 211; Table S1), indicate that the AlphaFold2-derived model (Fig. 1C) is resilient to large-scale changes due to mutagenesis. This stability is apparent when comparing the Ramachandran plots for the WT and variant simulations (Fig. S9A-F) [59]. The normalized B-factor values [58] indicate that most of the movement happen in the unstructured loops regardless of the system (Figs. S7 and S10). A re-adjustment can happen to accommodate the mutation; however, typically the RMSD curves level up at the end (Fig. S5).

A notable exception to this stability was the W572R variant, which underwent a major rearrangement in one of its replica simulations (Fig. 4A-C). Likewise, large-scale instability was induced with the WT system, when it was simulated at an elevated temperature (310 K → 400 K; Fig. S11). The truncated PH domain began to break its intra-protein interactions and rotate after ~20 ns in all three 400 K replica simulations. Because this ‘unfolding’ event could be repeated and it happened at the truncated domain, lacking the structural integrity to withstand the temperature boost, this was not studied further. Likewise, the simulation outcomes should not be considered reliable for the residues close to the model termini (‘Uncertain’; Table S1).

Although the missense mutations did not cause large-scale unfolding, they did cause notable conformational changes near the mutated residues (Figs. S5, S7, S9, S10). The residue swaps could affect the hydrophobic packing, disturb the intra-protein H-bonding network as well as weaken the secondary or tertiary structure. The variants were grouped into potentially benign (N = 52), potentially pathogenic (N = 117), and uncertain (N = 42) categories (Table S1). If the mutation was introduced near one of the truncated protein ends, the variant was given an uncertain status. No variant was assigned pathogenic based on its modelled location at the GAP-Ras or C2-membrane interfaces, unless the mutation had an adverse effect otherwise. A warning was given when there was a clear disparity between the expected impact and the simulation outcome (‘MD alert’; Table S1).

Mutations disrupting the secondary structure

Several missense mutations disrupt secondary structure elements in the SynGAP model based on the MD simulations (‘Secondary’; Table S1). Although these effects were often minor on the overall fold, due to the integral role of α helices and β sheets for the protein’s structural unity or the underlying folding process, the actual mutagenesis effects are potentially pathogenic [23, 60–64]. It is noteworthy that a secondary structure element can tolerate up to 20% of the ‘breaker’ residues (Gly, Pro, Ser, Asn, and Asp), and thus, a cluster of breakers can be needed to convert elements to loops or turns [65].

The negative effects on the secondary structure are straightforward with proline substitutions (Fig. 4D and Fig. S12A-B). Proline lacks the backbone amide needed for the H-bonding between the peptide bonds in the secondary structure elements [66–68]. In total, 14 proline variants, which clearly affect the integrity of secondary elements (e.g. R544P, L465P, and L431P), could be recognized. For example, in the L327P variant, the proline swap removes the backbone H-bond with the carbonyl of Gly344 at the middle of an anti-parallel β sheet in the C2 domain causing it to unfold (Fig. 4D). Likewise, with the L431P variant, the proline is introduced in the middle of an α helix, which removes the backbone H-bond with the carbonyl of His427 and leads to the helix breakage in the GAP domain (Fig. S12B).

Conversely, the rigid cyclic structure of proline facilitates tight peptide chain turns and limits the continuity of α helices and β sheets [69]. Although the proline replacements can affect the folding profoundly (e.g. P349S, P562L, P605R, and P605S), they did not generate drastic unfolding effects in the simulations. Similarly, glycine, second only to proline in low helix propensity [70], can transform an α helix into a loop due to its flexibility [71]; however, no such disruptions were visible (e.g. E419G, A541G, E496G, W572G, A577G, S590G, and E666G). Additionally, the serine swaps are known to generate bends or destabilize α helices, due to its sidechain hydroxyl H-bonding with the backbone atoms of the neighbouring residues [65, 66, 72]. Nevertheless, these ‘breakers’ did not affect the helix continuity (e.g. F420S, I510S, W572S, and L664S).

The β-hairpins, which consist of two antiparallel strands coupled by a turn and technically belong to super secondary structure motifs, have been proposed to act as nucleation sites for protein folding [73–78]. During annotation, missense mutations disrupting the β-hairpins were considered potentially pathogenic (e.g. L274Q, D287H, R299C, and L327P). The β-hairpins are abundant in the C2 domain that is composed of a β-sandwich (Fig. 1B). While in a few cases, the negative effect of the mutation on the hairpin was apparent (e.g. L327P; Fig. 4D), typically, it was difficult to ascertain their importance relying solely on the simulation outcomes.

Mutations affecting tertiary structure bonding

The SynGAP missense mutations that affect weak bonding interactions between domains or secondary structure elements, including flexible loops and turns, were considered as tertiary bonding effects (‘Tertiary bonds’; Table S1). The focus was given on gaining or losing H-bonds and salt bridges, because they can act as indicators of larger tertiary structure changes. Salt bridges and long-range electrostatic interactions are important for folding, facilitating tertiary assembly, and for improving the thermal stability [79–82]. The tertiary structure formation represents the final stages of the folding process [83], and, therefore, the variant simulations provide a viewpoint to the tertiary effects.

Indeed, several of the mutations were determined to affect the tertiary bonding or assembly based on the simulations (e.g. R342W, R405C, R429W, E436K, R573L, R596H, E666G, and R259W; Fig. S12C). For example, in the WT protein, Glu436 forms a salt bridge with Lys444 – a tertiary arrangement keeping two or even three α helices in contact. In the E436K variant, this inter-helix interaction is reversed by substituting the negative glutamate with a positive lysine that is repelled by the equally positive Lys444. This results in the shortening of the Lys444-containing helix and a compromised tertiary assembly (Fig. 4E). As a rule of thumb, when a charged group becomes either unpaired or wrongly paired via mutagenesis, as is the case with the E436K variant, misfolding can follow [65].

Inside-out mutations: Solvent exposure of buried residues

A missense mutation can render a buried residue unfit for its native intra-protein environment Table S1), which then leads to inside-out effects expulsing and/or misplacing sidechains. Because the hydrophobic effect is one of the major driving forces behind the protein folding [84], these inside-out effects can be detrimental to establishing a functioning protein. For example, with the C2 domain the inside-out mutations have a clear negative effect on the entire β-sandwich structure, although no specific bonds are necessarily altered near the mutation site (e.g. L264Q).

Mutations altering hydropathy, charge, and/or size, may expel the mutated sidechain (e.g. L323R, G390E, L402R, and A469D; Fig. 4F) or nearby residues from the protein interior (e.g. I494R, C547R; Fig. S12D-E). For example, the A496D variant introduces a negatively charged sidechain into the hydrophobic inter-helix space. The carboxylate of aspartate moves to the out-facing side of the α helix and forms a salt bridge with the guanidium of Arg575 to escape the electrostatically unfavourable environment, which, in turn, breaks the parent helix. Likewise, the I494R variant introduces a positive sidechain into a hydrophobic inter-helix space, where there is not enough room nor accommodating charge environment, which results in the loss of integrity of the α helix facing the introduced arginine.

Although the inside-out mutations are likely deleterious, typically no large-scale unfolding was seen in the simulations (e.g. W308R, C282R, and A271D; Fig. 5A and Fig. S12F, S12G). An exception to this was the W572R variant, where the expulsion of positive arginine from the inter-helix space caused the GAP and C2 domains to move apart (Fig. 4B and C); however, this unfolding event was not replicated.

SynGAP-Ras complex formation at the inner leaflet of the synaptic membrane

The WT-SynGAP-RasGTPase complex was modelled at the inner post-synaptic membrane based on four criteria to scrutinize its molecular interactions:

1. The C2 domain of SynGAP was determined to have Type II β sheet arrangement based on structural comparison against known Type I and II examples (Fig. S2B). The Type II C2 domains typically interact with the membrane surfaces in an upright position rather than via the sides, front or back of the β-sandwich fold [25]. Hence, the PH domain is not suggested to interact directly with the membrane, which is supported by the fact that most PH domains are not capable of independent membrane targeting [85].

2. The upright positioning of the SynGAP-C2 domain at the membrane surface was acquired by aligning it on top of the C2 domain of another protein (Fig. S2C) with settled membrane positioning [34, 35]. Notably, during the 150-ns WT-SynGAP-membrane simulations, the C2 domain remains stably membrane-bound, upright, and standing on ‘loop legs’ at the membrane surface (Fig. 2A).

3. The selected upright C2 domain positioning allowed several positive residues to face the polar or charged phospholipid head groups of the lipid bilayer in the SynGAP-membrane simulations (Fig. 2B and C; Table S2). Electrostatic forces have been suggested to be the main driving force behind the binding of Ca²⁺-independent C2 domains into membranes [25]. This not only supports the validity of the selected SynGAP-membrane assembly, but it also provides support to the notion that the binding of SynGAP on the inner membrane would follow the Ca²⁺-independent mechanism.

4. Lastly, the upright SynGAP-C2 domain positioning at the membrane is ideal for the GAP domain to bind unobstructed with its key interaction partner RasGTPase (Fig. 3A). The catalytic arginine finger of GAP binds to the Ras-bound substrate GTP (Fig. 3B) in the modelled SynGAP-Ras complex similarly as the equivalent domain of human GAP-Ras complex (Fig. 3C) [37]. The covalently linked Cys-palmitoylated lipids of Ras can easily dip into the inner membrane leaflet in this spatial arrangement. Ras has been postulated to have several different membrane orientations [86–91], it could even be adjusted by the nucleotide binding [88, 90], however, during dimerization [86] or when Ras is complexing with Raf [90, 91], a similar orientation has been suggested. Importantly, the WT-SynGAP-Ras complex was stably formed and membrane-bound in the simulation.

Mutations potentially altering the C2 domain-membrane association

From the perspective of the membrane-bound WT-SynGAP model (Fig. S2A), the C2 domain missense mutations – particularly those altering the membrane-facing loop dynamics (Table S1) – may compromise the stable SynGAP-membrane association. Here, the C2-membrane interaction was only probed with the WT protein; however, some interface mutations already displayed effects in the solvent simulations.

For example, with the A271D variant, the carboxylate of the aspartate is introduced into the C2 domain’s hydrophobic niche and, consequently, it forms H-bonds or salt bridges with the nearby loop residues (e.g. Lys394), which could decrease the flexible loop movement needed for the dynamic C2-membrane association (Fig. S12G). Likewise, with the V306D variant (Fig. 2D), the negative aspartate is repelled by the native hydrophobic partners of the WT residue valine (e.g. Trp308) and attracted by a new set of hydrophilic partners (e.g. Ser300). In general, the Ω-loops are known to have roles in protein functions mandating flexibility, and, thus, they rarely contain hydrophobic residues [92, 93]. Accordingly, the multiple valine substitutions into the Gly-rich Ω-loop motif of SynGAP (e.g. G373V, G381V, G387V, and G391V; Figs. 1B and 2D and Fig. S6) could reduce its flexibility, which, hypothetically, is needed for accomplishing the stable membrane association.

In theory, mutations that induce structural strain within the C2 domain (N = 85; Fig. S1) or disturb its direct interactions such as H-bonding with the membrane (Table S2; Fig. 2D), could compromise membrane binding and prove to be pathogenic. Nevertheless, in the absence of definitive evidence, these variants were only itemized, if the mutated residue was proximal to the membrane in the WT-SynGAP-membrane simulation.

Mutations potentially affecting the SynGAP-Ras association

The binding of the GAP domain to the RasGTPase (or RapGTPase) is a key effector interaction of SynGAP (Fig. 3A). Notably, its arginine finger (Fig. 3B) promotes the hydrolysis of GTP to GDP at the active site of Ras [61]. Although only the WT-SynGAP-Ras complex was simulated, from this perspective, any mutations directly at the protein–protein interface (Fig. 3D), for example, participating in the SynGAP-Ras H-bonding (Table S2; Fig. 3D) or relaying an allosteric effect from further away could perturb the control of GTPase activity.

There exists at least one missense mutation categorized as pathogenic (W362R; Table S1) and several VUS at the GAP-Ras interface that could weaken or even disrupt the protein–protein association (e.g. S604L, K642T, and R596C; N = 24; Table S1 and Table S2). For example, the hydrophobic sidechain of leucine in the S604L variant cannot H-bond with either Ras residue Ser65 or Ala66 similarly as the hydroxyl-containing serine of the WT protein. Likewise, in the K642T variant, the salt bridge formation between the WT residue lysine and the carboxylate groups of Asp33 and Asp38 at the Ras side is prevented by the residue substitution (Fig. 3D).

Because the connection between the compromised Ras activity and the SynGAP-related NSID is uncertain [6], variants that could alter the protein–protein interactions were only itemized (‘GAP-Ras interface’; Table S1) and given uncertain status (Table S1).

No structural effects: Potentially benign missense mutations

Plenty of SynGAP missense mutations have been categorized as harmless (N = 111; Table S1). While this benign or likely benign status can change due to emerging patient data, the MD simulations also suggest that many of the mutations are well-tolerated. Here, we considered a missense mutation potentially benign when it did not cause any apparent structural change, the effect was only minor, or its undesirable nature was unclear. For example, substituting valine with leucine (e.g. V400L; Fig. S12H), even when the sidechains are facing inward, the structural stability was not visibly affected during the simulations. A fraction of the missense mutations that were given the potentially benign status, are located at the protein surface, where they are unlikely to affect the fold or interact at least with the membrane or RasGTPase (e.g. A433V, V441A, L655Q, and W697R). Despite the uncertain status, several variants at the truncated model termini could be benign as well (e.g. D201E; Fig. 5A).

Pathogenic mutations destabilizing the structure

To complement the MD-based annotation, a quantitative modelling method was applied to explore the missense effects on the protein folding (Fig. 5B and C; Table S1). The FoldX-MD protocol [47], in which the mutations are assessed in the context of multiple WT simulation frames, is widely used to identify how single-point mutations affect the folding stability [47, 94–98]. The calculations predicted that 82 mutations destabilize the protein (e.g. C282R, W308R), 31 have a neutral effect (e.g. T532P, S535T), 89 fall into the ambiguous category, while none were predicted as stabilizing (Fig. 5C). The predictions should be considered unreliable for the PH region lacking structural integrity (Fig. S11). Although the calculations showed an 85% agreement with the MD-based annotation, a few differences were observed (e.g. D287Y, S300F, R429W, L465V, and G580S).

The folding stability calculations were useful in pointing out potentially pathogenic effects for the following mutation types:

1. The introduction of charged and bulky residues into the hydrophobic intra-protein positions was predicted to be strongly destabilizing (e.g. W308R; Fig. 5A).

2. Mutations that introduced a small polar residue like serine into the hydrophobic protein regions resulted in destabilizing energies (e.g. F420S, I510S, W572S, and L664S). Similarly, mutations that replaced a large hydrophobic residue with a smaller polar one, such as cysteine, were found to be destabilizing, in part due to the sidechain size difference (e.g. Y497C, W308C).

3. Variants that introduced proline within secondary structure elements were clearly unfavourable based on both the MD simulations and stability calculations (e.g. R293P, L323P, L431P, and L327P; Fig. 4D). However, the removal of proline from the ends of secondary structure elements via residue substitutions was also considered firmly destabilizing by the FoldX-MD calculations (e.g. P349S, P562L, P605R, and P605S) [67–69].

4. Residue swaps that introduced valine into the flexible Ω-loops were flagged as highly destabilizing (e.g. G373V, G381V, G387V, and G391V). The potential pathogenic effects of these mutations are understandable, provided that the C2-loops are in dynamic interaction with the membrane (Fig. 2A-C).

The FoldX-MD results were compared against variants with currently agreed upon non-VUS or non-contradictory status in the ClinVar to find out if the stability predictions matched the clinical data (Fig. 5B). Overall, the ΔΔG values indicate that the pathogenic (N = 11) or likely pathogenic variants (N = 31) were more often predicted to be destabilizing (e.g. W308R) than their benign (N = 17) or likely benign (N = 15) counterparts (e.g. D201E; Fig. 5A). Although this data set is limited (N = 74) and, by no means, should the status of the variants be considered settled, the comparison indicates an accuracy of ~ 87% Table S3) for the FoldX-MD calculations.

Because FoldX tends to predict stabilizing mutations with low accuracy [94, 99], Rosetta was also used to calculate folding energy [50, 99, 100]; however, it did not find any stabilizing mutations either (Table S1 available online at http://bib.oxfordjournals.org/). Rosetta had prediction accuracy of ~ 80% (Table S3 and Fig. S13A). Furthermore, Rosetta was more prone to assign benign status, while FoldX-MD was more effective at recognizing pathogenic mutations (Fig. 5A and Fig. S13A and Table S3). The methods showed ~ 73% overlap (Fig. S13D).