Identification of allosteric fingerprints of alpha-synuclein aggregates in matrix metalloprotease-1 and substrate-specific virtual screening with single molecule insights

Alpha-synuclein (aSyn) has implications in pathological protein aggregations in neurodegeneration. Matrix metalloproteases (MMPs) are broad-spectrum proteases and cleave aSyn, leading to aggregation. Previous reports showed that allosteric communications between the two domains of MMP1 on collagen fibril and fibrin depend on substrates, activity, and ligands. This paper reports quantification of allostery using single molecule measurements of MMP1 dynamics on aSyn-induced aggregates by calculating Forster Resonance Energy Transfer (FRET) between two dyes attached to the catalytic and hemopexin domains of MMP1. The two domains of MMP1 prefer open conformations that are inhibited by a single point mutation E219Q of MMP1 and tetracycline, an MMP inhibitor. A two-state Poisson process describes the interdomain dynamics, where the two states and kinetic rates of interconversion between them are obtained from histograms and autocorrelations of FRET values. Since a crystal structure of aSyn-bound MMP1 is unavailable, binding poses were predicted by molecular docking of MMP1 with aSyn using ClusPro. MMP1 dynamics were simulated using predicted binding poses and compared with the experimental interdomain dynamics to identify an appropriate pose. The selected aSyn-MMP1 binding pose near aSyn residue K45 was simulated and analyzed to define conformational changes at the catalytic site. Allosteric residues in aSyn-bound MMP1 exhibiting strong correlations with the catalytic motif residues were compared with allosteric residues in free MMP1, and aSyn-specific residues were identified. The allosteric residues in aSyn-bound MMP1 are K281, T283, G292, G327, L328, E329, R337, F343, G345, N346, Y348, G353, Q354, D363, Y365, S366, S367, F368, P371, R372, V374, K375, A379, F391, A394, R399, M414, F419, V426, and C466. Shannon entropy was defined to quantify MMP1 dynamics. Virtual screening was performed against a site on selected aSyn-MMP1 binding poses, which showed that lead molecules differ between free MMP1 and substrate-bound MMP1. Also, identifying aSyn-specific allosteric residues in MMP1 enabled further selection of lead molecules. In other words, virtual screening needs to take substrates into account for potential substrate-specific control of MMP1 activity in the future. Molecular understanding of interactions between MMP1 and aSyn-induced aggregates may open up the possibility of degrading aggregates by targeting MMPs.

Transmission Electron Microscope (TEM) images of aSyn-induced aggregates. We used freshly glow-discharged 200-mesh formvar carbon-coated copper grids (Electron Microscopy Sciences, Cat# FCF200-CU-50). We dipped two wooden tips in aSyn-induced aggregates and separated the two tips to create a thin layer of aSyn-induced aggregates on the TEM grid. After drying the sample for 5 min, we soaked the sample with 5 L of 2% uranyl acetate for 3 min and blotted any excess solution with a piece of Whatman filter paper. We rinsed the sample with 5 L of deionized water three times and let the sample dry for 5 min before TEM imaging. We imaged using a ThermoFisher Tecnai G2 Biotwin TEM (Cat# Tecnai12BT) at 80kV with an AMT sidemount XR80 8-megapixel digital camera. Figure S2 shows TEM images of aSyn-induced aggregates without MMP1 treatment. Note that we formed the aggregates and fibrils in the presence of other E. coli proteins, in contrast to fibrils formed with pure aSyn 1 . Best-fit parameters for experiments. We fitted a sum of two Gaussians to the experimental histograms in Figure 1. The best-fit parameters are in Table S1A. The fit parameters b1 and b2 are the two states, S1 and S2, respectively. We fitted exponential and power-law distributions to the experimental autocorrelations in Figure 1. Power law distribution does not fit the experimental autocorrelations. The best-fit parameters for the exponential fits are in Table S1B. The kinetic rates of interconversion between S1 and S2 are in Table S1C. The error bars represent the standard errors of the mean.

Calculation of correlations.
We used the following equation to calculate autocorrelations:  2, 3, 4, 1, 5, 3 3, 4, 1, 5, 3 We subtract the means from each element to obtain the fluctuations. We disregard the last element in the first row above because it is unpaired. For the first row, the mean is (2+3+4+1+5)/5=15/5=3. For the second row, the mean is (3+4+1+5+3)/5=16/5=3.2. After subtracting the means, we get the following time series: -1, 0, 1, -2, 2 -0.2, 0.8, -2.2, 1.8, -0. 2 We multiply element by element, add, and average by We need to shift the time series by the lag number 2 and multiply the time series element by element and calculate   C  at the lag number 2   as below: 2, 3, 4, 1, 5, 3 4, 1, 5, 3 We subtract the means from each element to obtain the fluctuations. We disregard the last two elements in the first row above because they are unpaired. For the first row, the mean is (2+3+4+1)/4=10/4=2.5. For the second row, the mean is (4+1+5+3)/4=13/4=3.25. After subtracting the means, we get the following time series: -0.5, 0.5, 1. We have followed this procedure to calculate normalized autocorrelations. For a two-state Poisson process, the decay of autocorrelations is the sum of two kinetic interconversion rates between the two states.
Best-fit parameters for simulations. We considered the MMP1 dynamics as a two-state Poisson process and simulated smFRET trajectories assuming that MMP1 undergoes interconversion between two states, S1 and S2, with experimentally determined kinetic rates and noise levels. We considered active MMP1 and active site mutant of MMP1 without ligands (Figures 1D and 1F). We simulated and analyzed 350 smFRET trajectories, each 1000 s long, with the input parameters in Table  S2. The recovered parameters are on the right side of Table S2. RMSD stabilization of simulations. We defined the input structure at t=0 as the reference structure and calculated the root-mean-square-displacement (RMSD) to check the simulations' stabilization.  Simulations of MMP1 dynamics stabilize in ~5 ns ( Figure S3). As such, we simulated 20 ns long dynamics for different conditions.

Calculation of entropy.
We calculated the Gray-Level Co-Occurrence Matrix (GLCM) 2 from the two-dimensional correlation plots. Then, we defined the Shannon entropy 3 from the GLCM. We describe the steps for calculating the Shannon entropy of MMP1 conformational dynamics in Figure S4 using an arbitrary 55 matrix, where each element of the matrix has a value between 0 and 3. In other words, we have a 2-bit 55 matrix in Figure S4. The dimension of GLCM depends on the number of possible values (0, 1, 2, 3), i.e., 44 for the 55 matrix in Figure S4. We calculated the GLCM at 0, i.e., we only considered the next neighbor on the left and right. For example, we calculated the (0,0) element of GLCM by finding the value 0s in the original 55 matrix and counting the number of times 0s appear on the left and right. The number is 0, and as such, the (0,0) element of GLCM is 0. To calculate the (0,1) element of GLCM, we found the value 0s in the original 55 matrix and counted the number of times 1s appear on the left and right. The number is 2, and as such, the (0,1) element of GLCM is 2. We repeated this process for all the elements of GLCM and obtained the matrix in the middle (Figure S4). We calculated the sum of all elements and divided each element by the sum to obtain the normalized GLCM at the right ( Figure S4) Figure S5. Metal Coordination Centre Parameterization. We used a bonded metal model approach to coordinate each metal ion with its corresponding protein coordination center. The six metal ions, two zinc and four calcium ions, with coordinating amino acids were modeled using the AmberTools21 MCPY.py (version 3.0) workflow. We advise the reader to consult the MCPB.py documentation for a complete step-by-step guide 4 . What follows is a brief overview of the steps and the parameters unique to the structure.

MD simulations with ions in
The complete coordination geometries of each metal ion coordination shell were first identified with a coordinating distance threshold of 2.8 angstroms using the 4AUO crystal structure of MMP1. Hydrogen atoms were added using the H++ web server 5 . The electronic structure of the extracted six coordination structures was optimized using GAMESS (US) 6 at a B3LYP/6-31g(d) level in a vacuum to remain consistent with earlier studies 7 . Charge and spin were set appropriately, and the protonation state of histidine was set according to its coordination with an ion. The normal modes and the frequencies were checked for imaginary numbers, and the Hessian was preserved in the output.
Force field parameters for Amber ff14SB were generated using the Seminario methods, and restrained electrostatic potential charges were fitted to each structure. The coordinates of a single MMP1 molecule from the 4AUO MMP1 crystal structure were used as a reference, and the residue names for all six metal ions and their respective coordinating side-chains were renamed to match the output files and parameter input file of MCPY.py. The Amber ParmEd program converted the final coordination and topology files suitable for the Gromacs suite. The parameters and topology data can be found at https://github.com/acnash/MMP1. Table S3. Best-fit parameters for simulated histograms and correlations in Figure 3.

Molecular Dynamics Parameters.
Atomic step integration was performed using the velocity Verlet algorithm with the time-step of 2 fs using Gromacs 2020.5 8 . The trajectory coordinates and velocities were recorded every 2500 steps (5 ps). Constraints were controlled using the LINCS schema, applied to hydrogen atoms only, with the LINCS iteration set to one with an order of four. The neighbor search was controlled using the Verlet cut-off scheme with the neighbor list updated every 10 steps. Periodic boundary conditions were applied in all directions. Long-range electrostatic interactions were calculated using the fast and smooth Particle-Mesh Ewald (PME) scheme, with a PME interpolation order of 4 and a Fourier spacing of 0.16 nm. A cut-off Van der Waals interactions and short-range Coulomb interactions were set to 1.0 nm, and both were modified with a potential shift. The temperature was set to 295 K for the two groups, water with counter ions and all protein and ligand-related atoms. The time constant for temperature coupling was set to 0.1 ps. Pressure coupling was set to 1 atmosphere and was controlled using the Berendsen followed by Parrinello-Rahman isotropic pressure coupling with a time constant of 2 ps. The compressibility was set to 4.5x10 -5 bar -1 . The reference coordinates were scaled using the center of mass. Table S4. Best-fit parameters for simulated histograms and autocorrelations in Figure 5.

Computational Model Construction.
Having performed force field parameterization of the six metal-ion coordination sites using the 4AUO crystal structure as explained, the experimental coordinates of each structure were encoding using the Amber ff14SB force field. Each structure was centered with a distance of at least 1.5 nm from every protein-ligand atom to the boundary of a dodecahedron unit cell, which was then solvated with tip3p water molecules with counter-ions for a neutral charge. Energy minimization using the steepest descent algorithm was performed with a stopping threshold of 800 kJ/mol. Position restraints of 1000 kJ/mol were applied to the alphacarbon atoms of the MMP-1 and alpha-synuclein ligand. Correlations are normalized between 0 (blue) and 1 (yellow). Yellowish colors indicate higher correlations. (F) Shannon entropy calculated from correlation plots for active (S= 3.05±0.01, mean±SEM) and inactive (S= 3.23±0.01, mean±SEM). For best-fit parameters, see Table S5.  Table S5. Best-fit parameters for simulated histograms and autocorrelations in Figure S4.
Each system was equilibrated for 50 ps using the NVT (fixed particles, volume, and temperature) ensemble at 295 K (22 C) using the velocity rescale scheme. The final frame coordinates and velocities were used for a further 50 ps while preserving the position restraints and switching to the Berendsen pressure coupling (NPT -fixed particles, pressure, and temperature -ensemble). Each simulation was continued for a further 50 ps with position restraints and a switch to the Parrinello-Rahman pressure coupling. Position restraints were removed, and each system was left to evolve for 100 ns. # do this many hybrid GA-LS runs analysis # perform a ranked cluster analysis Table S6. Best-fit parameters for correlation histograms in Figure 6.