We started from the X-ray crystal structure of the α2 I domain from α2β1 in complex with collagen [PDBid: 1dzi (Emsley et al., 2000)] and performed a scan of all possible mutations (for the 20 common amino acids) at each position along the “GFOGER” motif, collecting the expected free energy changes (ΔΔG) these programs predict. Integrin complexes were first optimized in the FoldX suite. Next, a position scan was conducted with the command “Position Scan” on the GFOGER motif and the output results showed the difference of binding energy for each mutation per amino acid on collagen. ΔΔG _bind was also calculated using RosettaDDG predictions, with the backrub trajectory stride set to 35,000 and making three trials for each ΔΔG calculation.

The ProteinMPNN (message passing neural network) (Dauparas et al., 2022) has recently been developed as a way to identify the ideal sequence that will adopt a certain 3D structure. In this model, we provided the PDB structure of the complex and asked the model to design new motifs to replace the native GFOGER motif.

We used standard minimization and equilibration protocols (Braun et al., 2018) followed by production runs using Langevin dynamics for 500 ns in the NPT ensemble using a Monte Carlo barostat (Åqvist et al., 2004). Simulations used AMBER’s (Case et al., 2020) pmemd module (Salomon-Ferrer et al., 2013). We simulated the top 20 FoldX and Rosetta predictions using ff14SB (Maier et al., 2015) solvated in a truncated octahedron box [OPC water model (Izadi et al., 2014)], and 150 mM concentration of Na⁺ and Cl⁻ ions (Joung and Cheatham, 2008). As a control, we simulated the I-domain in the presence and absence of the wild type (WT) collagen (PDBid: 1dzi). All simulations were carried out with a Co²⁺ ion in the MIDAS binding site. We simulated 10,000 steps of energy minimization, switching from steepest descent to conjugate gradient after 5,000 cycles. The resulting minimized system was heated from 0 to 100 K in NVT condition for 50 ps with Langevin dynamics, and 100–300 K in NPT for 500 ps using Langevin dynamics, followed by a short (5 ns) equilibration process at constant pressure (1 atm) and temperature (300 K). Finally, unbiased and unrestrained system went through production in a periodic boundary condition for 500 ns in NPT by Langevin thermostat and Monte Carlo barostat conditions. Bonds involving hydrogen were constrained by the SHAKE algorithm. Cpptraj (Roe and Cheatham, 2013) was used to analyze the root mean square deviation (RMSD) and Dynamical Cross Correlation (Kamberaj and Vaart, 2009) within the ensembles comparing them to the wild type complex.

We used Alphafold Multimer (Evans et al., 2021b) to predict the structure of the complex using either sequence data or templates (containing the collagen and integrin domain far from each other). Results were analyzed in terms of the predicted local distance difference test (pLDDT) score as is standard in the field (Jumper et al., 2021). In short, the pLDDT score gives a per residue and global value to show how confident the Alpha Fold prediction results are. Results above 80 typically reflect high confidence in the prediction.

TI was used to calculate the relative binding affinity ( Δ Δ G bind mutant−WT ) between collagen and α2β1 upon mutation of certain residues in collagen. Here, we applied “One-step” transformations (Steinbrecher et al., 2011) to decrease the simulation time, in which electrostatic and van der Waals forces are varied synchronously (Shirts et al., 2003). The initial system was prepared using AMBER’s tiMerge to eliminate redundant bonding terms and increase calculation efficiency. We ran TI simulations with pmemd. The complex and mutant ligand were solvated separately in a cubic box with explicit OPC (Izadi et al., 2014) water and a 10 Å clearance. We employed ff14SB (Maier et al., 2015) for the protein parameters and general AMBER force field (He et al., 2020) for general atom and bonds parameters. Minimization, heating and equilibrium process was performed in the NVT ensemble with a Monte Carlo barostat. The TI production phase was done in the NPT ensemble (300 K and 1 atm), running for 500 ns Softcore potentials were applied to reduce issues with the integration step at the endpoints (Steinbrecher et al., 2011). Eleven independent MD simulations were performed spaced evenly between the end-points (λ ∈ [0, 1]). We performed six replicates for each simulated system. The average and standard deviation for ΔΔG _bind were calculated from the differences amongst replicates.

The Modeling Employing Limited Data (MELD) approach uses H,T-REMD (Sugita and Okamoto, 1999) to sample rare events. The method changes the Hamiltonian by enforcing information that guides to different conformations that might be compatible with the end state. The caveat is that the data is framed as ambiguous and noisy—thus MELD relies on Bayesian inference to identify the best interpretation of the data compatible with the forcefield. In this process, analyzing the resulting ensemble (e.g., through clustering) identifies the states (conformations) most compatible with the information and force field.

To guide the binding process we first placed harmonic distance restraints amongst native contacts in the integrin (so it would not unfold), and also between the three collagen strands, so it would not dissociate. We then selected residues in the active site of the integrin and in those of the collagen binding motif. Based on those two lists of residues, we generated a list of twenty five possible contacts (some of which were present in the native state and some of which were not). We found that when enforcing 15 or more restraints, replica exchanges were inefficient, leading to poor sampling. At the other extreme, satisfying less than four restraints sampling was not restrictive enough to sample native-like bound conformations. We thus required that only eight restraints out of the 25 possible ones be satisfied. Satisfying different subsets of eight restraints give rise to different binding modes.

MELD simulations used the ff14SB force field (Maier et al., 2015) for side chains and ff99SB (Hornak et al., 2006) for backbone, together with the GBneck2 (Nguyen et al., 2013) implicit solvent model. The collagen fiber was placed over 30 Å away from the integrin. The temperature range was set between 300 and 500 K, with 30 replicas. Ensembles were analyzed using hierarchical clustering as implemented in CPPTRAJ (Roe and Cheatham, 2013) with an ϵ = 2 value, including heavy atoms at the interface of the complex in the native state.

Traditional design strategies start with a known binding motif and search for single amino acid mutants that increase binding affinity (ΔΔG _bind ). Such strategies lead to local sequence optimization, with designs similar to the original motif. Here we used FoldX (Schymkowitz et al., 2005) and Rosetta (Barlow et al., 2018) (see methods), two traditional approaches with varying computational cost and success rate. We observed that FoldX single point mutations have a wider ΔΔG _bind distribution, and are generally shifted towards higher energies (see Figure 2A). While there is good agreement on the failed mutations, the more computationally demanding Rosetta is better at discriminating mutations that FoldX finds favorable.

To further assess the predicted motifs with an independent methodology, we performed MD simulations of a selected group of 15 mutants. We expect that monitoring standard structural and dynamical properties like RMSD and dynamical cross correlation functions would be enough to distinguish those mutations that remain stable in the 500 ns timescale vs. those that are unlikely to bind (see Figures 2C–F). We monitored the RMSD of the interface region, defined as heavy atom contacts to collagen in the native structure (using a 10 Å cutoff). In this timescale, the integrin oscillates around 1 Å from the initial structures, with few deviations to higher RMSD values (2.5 Å). In the presence of collagen we observe a similar behavior, where there are no deviations to larger RMSD states in the 500 ns timescale. The RMSD of the whole complex oscillates at around 4 Å. Figure 2 showcases the behavior of the wild type, neutral and negative mutation on sampling [RMSD and projection on the top two principal components using the Bio3d package (Grant et al., 2006)]. Figure 2 exemplifies a negative control mutation (which rapidly dissociates) and a neutral mutation that remains close to the starting conformation.

We find that the more computationally efficient FoldX is capable of filtering out mutations that are likely detrimental to the binding affinity. While the ones predicted to be beneficial do not always agree with MD and Rosetta results (see Figure 2B). We notice several disagreements with Rosetta and MD—this is not surprising as Rosetta has been designed to predict free energy differences while short conventional MD trajectories do not contain enough sampling to assess the free energy. We thus decided to perform thermodynamic integration calculations to further identify the agreement between Rosetta and MD-based approaches.

Thermodynamic integration increases the complexity in system setup and analysis with respect a conventional MD trajectory—but the computational costs (considering replicates needed, see methods) remains relatively small compared to other MD approaches. We selected 15 mutations and compared results using Rosetta and TI (see Supplementary Figure S1). For most residue mutations, both programs agree in sign if not in magnitude. Previous work points to systems including multiple binding modes or systems that are sensitive to local conformational changes (such as the MIDAS binding site) (Armacost et al., 2020) as problematic for TI. For example, Guest and coworkers performed free energy perturbation studies on a series of small molecule inhibitors to the β6 integrin with an average error of 1.5 kcal/mol with respect to the experimental results (Guest et al., 2020).

We searched for alternative binding modes by using the MELD approach, which can simulate multiple binding/unbinding events. MELD combines ambiguous/noisy information with molecular simulations through Bayesian inference and has been routinely used for predict the binding of macromolecules [protein-protein (Brini et al., 2019), protein-peptide (Morrone et al., 2017; Mondal et al., 2022), protein-DNA (Bauzá and Pérez, 2021), and protein-small molecule (Liu et al., 2020)]. We derived ambiguous information based on native contacts present in the crystal structure in such a way that different interpretations of the data is compatible with different binding modes. We expected, that the force field would be able to recognize the most native-like amongst the binding modes for those sequences that have a high affinity (clusters with high population) (Lang and Perez, 2021). Unfortunately, due to the small interface region between collagen and the integrin, the different binding modes found give rise to large deviations in binding angles between the collagen in MELD simulations with respect to the native structure (see Supplementary Figure S2). On the other hand, satisfying more information overrides the force field preferences and yields native-like binding modes regardless of the sequence. Similarly, competitive binding simulations (Morrone et al., 2017) with MELD also failed to distinguish which collagen mutations were more likely to lead to more stable complexes. Presumably, these limitations arise from the use of an implicit solvent (Nguyen et al., 2013) needed for the MELD binding simulations.

Similarly, the recent successes of the AlphaFold (AF) (Evans et al., 2021a) machine learning approach did not translate to this system. We used a local installation of AlphaFold and performed predictions in the presence/absence of structural templates. In our hands, Alphafold multimer predictions were confident about the α2 I-domain structure (high pLDDT scores), but failed to predict the structure of the collagen triple helix structure—and hence of the complex (see Supplementary Figure S2).

Whereas we used FoldX and Rosetta to predict local changes in the sequence (single mutants), the recent protein MPNN (Dauparas et al., 2022) machine learning approach can in principle find an optimal sequence given the structure of the complex. Contrary to the other two methods, this approach does not provide a relative binding affinity. We first generated two predictions in which we allowed any residue in the motif along the tree collagen strands to change (see “Prediction 1″ and “Prediction 2” in Supplementary Figure S3). This gave rise to four different binding motifs. We next generated four more sequences by creating homo-trimer collagen strands with each of the four predicted motifs (see the latter four motifs in Supplementary Figure S3A). We assessed the viability of these motifs by running conventional MD. All sequences in which the leading strand had an E to P mutation were unstable. Whereas if this mutation occurred in other strands, the system remained stable. This is expected as the Glutamic acid coordinates with a divalent site when interacting with the integrin.