To probe the mechanics of the interaction between Bbp and Fgα, and to characterize the atomic details of the complex under force load, we performed aa-SMD simulations with Bbp anchored by its C-terminal while Fgα was pulled at different velocities (Supplementay Table S1). The simulations resulted in Force vs extension curves that reveal a clear one-step rupture event, as represented in Figure 2A. For the slowest pulling velocity, 160 replicas were performed following a wide-sampling paradigm previously developed in our group Sedlak et al. (2020). At the pulling velocity of 2.5 × 10⁻⁰⁴ nm/ps, we observed that the most probable rupture force for the complex was 3,510 pN, as described by the Bell-Evans (BE) Bell (1978); Evans and Ritchie (1997) fit of the peak forces at that pulling speed (see Figure 2B). Our results reveal that Bbp: Fgα is the most mechanostable complex investigated thus far, which is in agreement with previous experimental data where we showed that SdrG: Fgβ complex can withstand forces on the 2 nN range, equivalent to breaking of covalent bonds Milles et al. (2018).

To investigate the dependence of the mechanostability of Bbp: Fgα on the force loading rate, we performed CG-SMD simulations at several, much lower, pulling speeds (Supplementary Table S1). We have recently shown that aa-SMD and CG-SMD can be combined to in an in silico SMFS approach Gomes D. E. et al. (2022); Melo et al. (2022). Here, the combination of the two levels of molecular details is capable of rendering predictions that are consistent with theory and experimentation with the advantage of being 10 to approximately 100 times faster than aa-SMD simulations, depending on the pulling speed Gomes D. E. et al. (2022); Melo et al. (2022). A Dudko-Hummer-Szabo Dudko et al. (2006) (DHS) fit was performed through the SMD data, including both the aa-SMD, and the CG-SMD (see Figure 3). The DHS fit suggests that the system should rupture at forces higher than 2 nN at 10⁵ pN/s force loading rate, in agreement with experimental data Milles et al. (2018). It is interesting to note that the BE model is able to fit well all the simulation results, at both aa and CG level, as evidenced by the density plots in Figure 3.

The influence of the peptide size on the rupture force was also investigated. We have shown previously that SdrG complexed with shortened Fgβ peptides had lower unbinding forces Milles et al. (2018). Here, we simulated a model of Bbp complexed with Fgα elongated by nine residues (See Methods section) by aa and CG-SMD simulations (Supplementary Table S1). Our results show that the force loading rate was not significantly impacted by the size of the peptide (Supplementary Figure S1), indicating that the original complex formed at the crystal structure has the minimal length to keep the important contacts with the protein latch to hold the DLL configuration.

After confirming that Bbp: Fgα complex presents a hyperstable interaction under shear mechanical load, we used the approximately 3 μs of aa-SMD simulation data to investigate the molecular origin of the mechanostability of the complex. Previously, simulations of the SdrG: Fgβ revealed the presence of frequent and persistent hydrogen bonds (H-bonds) between the peptide and the protein backbone, showing that the high-force resilience of the complex was largely independent of the peptide side-chains interactions, and therefore the peptide’s sequence Milles et al. (2018). Here, we computed the occupancy of the H-bonds between the Bbp and Fgα before the complex rupture. We identified the key amino acid interactions responsible for keeping the complex together at high force loads (Table 1). Different than SdrG: Fgβ, Bbp: Fgα interactions are not dominated by backbone-backbone interactions, with a significant amount of side-chain interaction of the peptide playing an important role in the complex mechanostability. The backbone interactions between Bbp ^{Leu584, Thr582, Thr586} and Fgα ^{Thr565,Ser567,Thr586} have been previously described as important for Fgα binding at the crystal structure Zhang et al. (2015). However, we noticed that the side-chain H-bonds are rearranged upon application of mechanical stress on the complex. On the crystal, Bbp^Asp334 forms a side-chain H-bond with Fgα ^Ser566, and during the SMD simulations, this interaction shifts to Fgα ^Thr565, being the H-bond with the highest occupancy over the trajectories. Another shift occurs between Bbp^{Asp334,Ile335} interacting with Fgα ^Phe564, on the crystal, to Bbp^Ser333 interacting with Fgα ^Phe564 in our simulations. The H-bond between Bbp^Asp588 and Fgα ^Gln563 is described as important to lock the peptide N-terminus and is still present before the rupture of the complex, although with lower occupancy. Instead, a charged side-chain interaction arises with significant occupancy values: Bbp^Asp556:Fgα ^Lys562. These data corroborates the importance of backbone interactions to maintain the high mechanostability and also highlights important side chain H-bonds plasticity that occurs when Bbp: Fgα is exposed to mechanical stress.

How a shear force load “activates” the hyperstability of the complex can be investigated by analysing the evolution of pairwise interactions during the force-loading event. Such analysis can be used to investigate how a catch-bond may be formed in the Bbp: Fgα complex Liu et al. (2020). Previously, it has been shown that SdrG: Fgβ presents a catch-bond behavior Huang et al. (2022), which is expected also for Bbp: Fgα. To analyse the pairwise interactions during the SMD, we employed the generalized correlation-based dynamical network analysis method Melo et al. (2020), which can also be used to calculate force propagation pathways Schoeler et al. (2015). Figure 4A shows the pairwise interactions obtained from the network analysis. The thickness of the connections between nodes (amino acid residues) represents how well correlated the motion of these nodes are, and therefore how well connected are these amino acid residues.

The force propagation pathway that connects the pulling and the anchoring residues shows that most of the force is propagating from the protein latch directly to the peptide, passing by the center of Bbp’s N2 domain (Figure 4B). These results are slightly different than the ones obtained for the SdrG: Fgβ complex upon high mechanical stress Milles et al. (2018). However, in a previous study, it was observed that changes in the pulling velocities can lead to different force propagation pathways, suggesting different unbinding mechanisms at different pulling rates Melo et al. (2022).

The rigidity of the protein under high-force load can also be studied using the betweenness map from the dynamical network analysis (see Figure 4C). The betweenness is defined as the number of shortest paths from all vertices to all others that pass through that node, in this case, an amino acid residue. If an amino acid residue has high betweenness, it tends to be important for controlling inter-domain communication within a protein Melo et al. (2020). High betweenness values (thicker red tubes) are seen on the latch that is in direct contact with Fgα, highlighting the strong correlation between the motif and the peptide. Interestingly, high betweenness is also found at connections intra N2 domain, pointing that Bbp: Fgα complex becomes more rigid under high force loads, particularly in the region interconnecting the latch, the peptide and the N2 domains. Such behavior helps the stabilization of the interactions under high forces.

A representation of the network in subgroups, or communities, is shown at Figure 4D. The communities group the amino acid residues that are most inter-connected in relation to the rest of the network. We can see that Bbp: Fgα is subdivided in a handful of communities. The latch, most of Fgα, and part of the N2 domains are united in the same community in light blue, showing that these amino acids are highly connected. We also measured the correlation between motions on the interface residues to determine how cooperative their motion is and the essential contacts that are keeping the complex stable under high mechanical load. Essentially, the higher the correlation between residues, the more relevant is their interaction for the stability of the protein complex. We noticed that two Fgα residues are highly correlated (values equal or superior to 0.5) to Bbp at the interface, namely: Fgα ^Gln563: Bbp^{Asp588,Ser585,Thr586,Thr587} and Fgα ^Phe564:Bbp^Ser585 (Figure 5). The importance of Fgα ^Gln563 described as a persistent H-bond contact with Bbp^{Asp588,Thr586} and important locking contact is reinforced by its high correlation values. The same analysis was performed for the simulations of Bbp complexed with the elongated Fgα peptide (Supplementary Figure S2). The pattern of contacts is very similar, reinforcing the importance of Fgα ^Gln563, and we observe the absence of new contacts made by the extra residues, corroborating that the short peptide contains the key residues responsible for holding the complex tight at the DLL configuration.

The structure of Bbp in complex with Fgα has been previously solved by means of X-ray crystallography at 1.45 Å resolution Zhang et al. (2015) and deposited at the Protein Data Bank (PDB ID: 5CFA). Here we retrieved this structure and prepared it for molecular dynamics (MD) simulations using VMD Humphrey et al. (1996) and its plugin QwikMD Ribeiro et al. (2016). To investigate the loading rate dependency on the size of the peptide, we used Modeller v.10.1 Webb and Sali (2016) to create an additional structure of the complex where the Fgα was elongated by nine residues at its C-terminal end, in respect of the crystal structure, following the sequence of Fgα from Homo sapiens (Uniprot ID: P02671). The model followed the same preparation as described for the crystal structure.

The complexes between BBP and Fgα in its short or longer configuration were solvated using the TIP3P water model Jorgensen et al. (1983), with the net charge of the protein neutralized using a 150 mM concentration of sodium chloride. Steered molecular dynamics (SMD) simulations were carried out using NAMD 3 Phillips et al. (2020), with the CHARMM36 force field Best et al. (2012). The simulations were performed assuming periodic boundary conditions in the isothermal-isobaric ensemble (NPT) with temperature maintained at 300 K using Langevin dynamics for temperature and pressure coupling, the latter kept at 1 bar. A distance cut-off of 11.0 Å was applied to short-range non-bonded interactions, whereas long-range electrostatic interactions were treated using the particle-mesh Ewald (PME) Darden et al. (1993) method. Taking advantage of a hydrogen-mass repartitioning method implemented in VMD’s autopsfgen, the time step of integration was chosen to be 4 fs for all production aa-MD simulations performed. Before the SMD simulations, the system was submitted to an energy minimization protocol for 1,000 steps. An MD simulation with position restraints in the protein backbone atoms was performed for 1 ns, with temperature ramping from 0 K to 300 K in the first 0.5 ns at a timestep of 2.0 fs in the NVT ensemble, which served to pre-equilibrate the system. In an in silico single-molecule force spectroscopy (SMFS) strategy Verdorfer et al. (2017); Bernardi et al. (2019), SMD simulations were carried out in several replicas, using a constant velocity stretching protocol at three different pulling speeds (Supplementary Table S1). SMD was employed by harmonically restraining the position of the amino acid at the C-ter of Bbp and moving a second restraint point at the C-ter of Fgα peptide with a 5 kcal/mol Ų spring constant, with constant velocity in the z-axis. The force applied to the harmonic spring is then monitored during the time of the SMD. The pulling point was moved with constant velocity along the z-axis and due to the single anchoring point and the single pulling point the system is quickly aligned along the z-axis. The number of replicas for each velocity is indicated at Supplementary Table S1.

The atomistic model of Bbp: Fgα was modeled onto the Martini 3.0 Coarse-grained (CG) force field (v.3.0.b.3.2) Souza et al. (2021) using martinize2 v0.7.3 Kroon (2020). A set of native contacts, based on the rCSU + OV contact map protocol, was computed from the equilibrated all-atom structure using the rCSU server Wołek et al. (2015) and used to determine Gö-MARTINI interactions Poma et al. (2017) used to restraint the secondary and tertiary structures with the effective depth ϵ of Lennard-Jones potential set to 9.414 kJ.mol⁻¹. All CG-MD simulations were performed using GROMACS version 2021.5 Abraham et al. (2015). The Bbp: Fgα complex was centered in a rectangular box measuring with 10.0, 10.0, 25.0 nm to the x, y, and z directions. The anchor (Bbp C-terminal) and pulling (peptide C-terminal) backbone (BB) atoms were used to align the protein to the Z-axis. The box was then solvated with Martini3 water molecules. Systems were minimized for 10,000 steps with steepest descent, followed by a 10 ns equilibration at the NPT ensemble using the Berendsen thermostat at 298K, while pressure was kept at 1 bar with compressibility set to 3e⁻⁴bar⁻¹, using the Berendsen barostat. A time step of 10 fs was used to integrate the equations of motion. Pulling simulations were subsequently done at the NVT ensemble with a time step of 20 fs the temperature was controlled using the v-rescale thermostat Bussi et al. (2007) with a coupling time of 1 ps for all CG-MD simulations, the cutoff distance for Coulombic and Lennard-Jones interactions was set to 1.1 nm De Jong et al. (2016), with the long-range Coulomb interactions treated by a reaction field (RF) Tironi et al. (1995) with ϵ _r = 15. The Verlet neighbor search Verlet (1967) was used in combination with the neighbor list, updated every 20 steps. The LINCS Hess et al. (1997) algorithm was used to constrain the bonds and the leapfrog integration algorithm for the solution of the equations of motion. Several replicas of CG-SMD simulations were performed at a range of speeds described at Supplementary Table S1.

All analysis presented at the main text correspond to the Bbp: Fgα original complex. Force loading rate and mean correlation values for Bbp complexed with the elongated Fgα peptide are found on the Supplementary Information material. H-bonds occupancy between Bbp and Fgα were calculated and averaged for aa-MD simulations 1 ns before the main rupture event, using VMD Humphrey et al. (1996) with standard parameters for the calculation: residue pairs; donor-acceptor distance of 3.0 Å; angle cutoff of 20°. Mean correlation and dynamical network pathways were calculated using the generalized dynamical network analysis Melo et al. (2020) and VMD for aa-SMD at pulling velocity of 2.5 × 10⁻⁴ nm/ps. In this analysis, a network is defined as a set of nodes that represent amino acid residues, and the node’s position is mapped to that of the residue’s α-carbon. Edges connect pairs of nodes if their corresponding residues are in contact, and two non-consecutive residues are said to be in contact if they are within 4.5 Å of each other for at least 75% of analyzed frames. The interface residues between Bbp: Fgα were defined in a radius of 10 Å between nodes in each molecule. A representative for the full-network, optimal and suboptimal paths and communities was rendered using one of the SMD trajectory replicas. The mean correlation analysis was carried out 1 ns before the first rupture event using a cutoff of 0.35 for the mean correlation coefficients. All charts were generated using in-house python scripts. The protein image was rendered using VMD.