4.5. SILCS Simulation

Site Identification by Ligand Competitive Saturation (SILCS) simulation is used to explore functional group affinity patterns on both Ca²⁺-S100B and Ca²⁺-S100A1. The SILCS method involves molecular simulations of the target protein immersed in an aqueous solution that contains additional organic solutes of different chemical classes [37]. The solutes and water then compete for binding sites on the protein surface during the simulation, yielding a free energy fragment competition assay from which 3D fragment probability distributions of the solutes are used to define affinity patterns, termed FragMaps, encompassing a dynamic protein surface.

The current SILCS run was performed using the oscillating chemical potential Grand Canonical Monte Carlo (GCMC)/MD protocol for SILCS [38,39]. The structure models of S100B and S100A1 developed as described above was used to initialize the SILCS simulations. The protein was solvated in a water box, the size of which was determined to have the protein extrema separated from the box edge by 12 Å on all sides. Eight representative solutes with different chemical functionalities (benzene, propane, acetaldehyde, methanol, formamide, imidazole, acetate, and methylammonium) were added into the system at ~0.25 M concentration, to probe the functional group requirements of the protein. Ten such systems with different fragment positions were created to expedite the convergence of the simulations. Each system was minimized for 5000 steps with the steepest descent (SD) algorithm in the presence of periodic boundary conditions (PBC) and was followed by a 250 ps MD equilibration [40]. During SILCS simulations, weak restraints were applied on the backbone alpha carbon atoms with a force constant (k in 1/2 kδx²) of 0.12 kcal/mol/Ų to limit large conformational changes in the protein and to prevent the rotation of the protein in the simulation box. The ten GCMC/MD simulations were run for 100 cycles where each cycle has 200,000 steps of GCMC and 1.0 ns of MD, yielding a cumulative 200 million steps of GCMC and 500 ns of MD. During GCMC, solutes and water are exchanged between their gas-phase reservoirs; the excess chemical potential used to drive such exchange is varied every 3 cycles to yield an average concentration corresponding to 0.25 M of each fragment [38]. The configuration at the end of each GCMC run is used as the starting configuration for the following MD. During MD, the Nosé−Hoover method was used to maintain the temperature at 298 K and pressure was maintained at 1 bar using the Parrinello−Rahman barostat [41,42,43]. CHARMM36m protein force field, CHARMM General Force Field (CGenFF) and modified TIP3P water model [44,45,46,47] were used to describe protein, organic solutes, and water during the simulation, respectively. GCMC was performed and MD was conducted using GROMACS program [39,48].

3D probability distributions of the selected atoms from the organic solutes, called “FragMaps,” from the SILCS simulations were constructed and combined to obtain both specific and generic FragMap types as previously described [49]. Atoms from snapshots output every 10 ps from each SILCS simulation trajectory were binned into 1 Å × 1 Å × 1 Å cubic volume elements (voxels) of a grid spanning the entire system to acquire the voxel occupancy for each FragMap atom type being counted. The voxel occupancies computed in the presence of the protein were divided by the value in bulk to obtain a normalized probability. Normalized distributions were then converted to grid free energy (GFE) based on a Boltzmann transformation for quantitative use [49].