where KD is the apparent binding dissociation constant, [L]T and [P]T are the total ligand and protein concentrations used in the experiment and ∆δobs is the observed change in chemical shift upon ligand addition. ∆δu and ∆δb represent the changes in chemical shift for the unbound and bound states, respectively. Based on the KD values and accounting for fitting error, residues that were interacting with the ligands were identified and clustered into binding sites. Curve fitting and calculations were performed on Matlab R2019b and the protein surface visualization was carried out on PyMol 2.3.5 viewer (Schrödinger).

Protein Surface Properties

A structural model for the IgG1 FC domain was built by homology modeling starting with the crystal structure for an aglycosylated human IgG1 FC fragment (PDB code: 3S7G) using Molecular Operating Environment (MOE 2018, Chemical Computing Group). The electrostatic potential (EP) maps were calculated using the adaptive Poisson-Boltzmann solver (APBS) (Baker, Sept, Joseph, Holst, & Andrew McCammon, 2001) and surface aggregation propensity (SAP) maps were calculated as described by Trout and co-workers (Chennamsetty, Voynov, Kayser, Helk, & Trout, 2010). The EP map was calculated at pH 5.0 to match the experimental conditions. The resulting protein surface maps were then visualized using the PyMol 2.0.6 viewer (Schrödinger).

Molecular Dynamics Simulations

MD simulations were performed with the two MM CEX ligands, Capto MMC and Nuvia cPrime in free solution around the FC molecule, as shown in Figure 2. Each simulation was performed for 200ns with a timestep of 2fs and storing one frame every 1ps. The first 50ns were taken as equilibration time and all analyses were performed using the last 150ns. The simulation box dimensions were 8.5nm x 10nm x 13nm, allowing for a buffer of roughly 1.5nm on each side of the protein. The FC was prevented from rotating in each simulation by restraining a single alpha carbon buried in the center of each of the four Fc domains using a harmonic potential with a spring constant of 40,000kJ/mol/nm2 (Srinivasan et al., 2017). This allowed the use of a rectangular simulation box without risking the protein interacting with itself through periodic boundary conditions. Additionally, each simulation contained sodium counterions for electroneutrality and an excess of 18 sodium and chloride such that there is a counterion for every charged side chain. A total of 61 ligands were included in each simulation, corresponding to a concentration of 0.1 M. The FC was parameterized using the AMBER Parm99 (Cheatham, Cieplak, & Kollman, 1999) forcefield and PropKa (Bas, Rogers, & Jensen, 2008) was used to adjust the charge of the side chains to reflect a pH of 5.0. Water was modeled explicitly using the TIP3P (Jorgensen et al., 1981) water model and ligand atom types and nonbonded interactions were parameterized using the General AMBER Force Field (GAFF) (Wang, Wolf, Caldwell, Kollman, & Case, 2004). Atomic partial charges were obtained using a GAUSSIAN (Frisch, M.J.E.A., Trucks, G.W., Schlegel, H.B., Scuseria, G.E., Robb, M.A., Cheeseman, J.R., Scalmani, G., Barone, V., Mennucci, B., Petersson, 2009) calculation with RESP (Bayly, Cieplak, Cornell, & Kollman, 1993; Wang, Cieplak, & Kollman, 2000) assignment using the Antechamber tool of AMBER (Wang, Wang, Kollman, & Case, 2001), as has been described previously (Bilodeau, Lau, Cramer, & Garde, 2019). The resulting topologies were converted into the GROMACS format using ACPYPE (Silva & Vranken, 2012). In this work, partial negative charges on the carbons in the phenyl ring were employed rather than explicit parameterizations to reflect pi interactions.
Simulations were performed using GROMACS 4.5.3 (Hess, Kutzner, Van Der Spoel, & Lindahl, 2008; Pronk et al., 2013) in the NPT ensemble. The temperature and pressure were maintained at 298 K and 1 atm using a Nosé-Hoover thermostat (Evans & Holian, 1985) and a Parrinello-Rahman barostat (Parinello, M. and Rahman, 1980), respectively. Ligands, protein, and water/ions were treated as three separate temperature coupling groups. Electrostatic interactions were calculated using the Particle-Mesh Ewald (Darden, York, & Pedersen, 1993) method with a grid spacing of 0.1 nm, an order of 4 for the B-spline interpolation, and a direct sum tolerance of 10-5 (consistent with default parameters).