Application of Funnel Metadynamics to the Platelet Integrin αIIbβ3 in Complex with an RGD Peptide

Integrin αIIbβ3 mediates platelet aggregation by binding the Arginyl-Glycyl-Aspartic acid (RGD) sequence of fibrinogen. RGD binding occurs at a site topographically proximal to the αIIb and β3 subunits, promoting the conformational activation of the receptor from bent to extended states. While several experimental approaches have characterized RGD binding to αIIbβ3 integrin, applying computational methods has been significantly more challenging due to limited sampling and the need for a priori information regarding the interactions between the RGD peptide and integrin. In this study, we employed all-atom simulations using funnel metadynamics (FM) to evaluate the interactions of an RGD peptide with the αIIb and β3 subunits of integrin. FM incorporates an external history-dependent potential on selected degrees of freedom while applying a funnel-shaped restraint potential to limit RGD exploration of the unbound state. Furthermore, it does not require a priori information about the interactions, enhancing the sampling at a low computational cost. Our FM simulations reveal significant molecular changes in the β3 subunit of integrin upon RGD binding and provide a free-energy landscape with a low-energy binding mode surrounded by higher-energy prebinding states. The strong agreement between previous experimental and computational data and our results highlights the reliability of FM as a method for studying dynamic interactions of complex systems such as integrin.


Introduction
Integrins are transmembrane heterodimers that play a crucial role in mediating adhesions and signaling pathways between cells and their surrounding environment or other cells.Comprising α and β subunits (Figure 1A,B), there are currently 18 α and 8 β subunits known, forming a total of 24 different combinations [1].Among these, the α IIb β 3 heterodimer is expressed by blood platelets and is essential for hemostasis, platelet aggregation, and blood clot retraction, making it a significant target for antithrombotic therapies.One of the best-established paradigms in α IIb β 3 integrin biology is the recognition of ligands through the Arg-Gly-Asp (RGD) sequence (Figure 1C,D).The affinity of α IIb β 3 for the RGD sequence increases as conformational activation progresses from bent-closed to extended-open states [2,3].Additionally, α IIb β 3 recognizes an AGDV sequence at the C-terminus of the fibrinogen γ subunit [4], and once fibrinogen is converted to fibrin by thrombin, integrin binds other sites composed of RGD sequences [5][6][7].RGD mimetic drugs compete with fibrinogen binding to α IIb β 3 .Currently, three RGD-based drugs, namely tirofiban, abciximab [8], and eptifibatide [9], are clinically used to inhibit platelet aggregation and prevent thrombosis.Their crystal structures complexed with α IIb β 3 have been determined [10].However, these drugs act as partial agonists [11,12], leading to increased thrombosis in some patients and exposure to neo-epitopes, causing an immune reaction that results in thrombocytopenia in others [13].RGD ligand mimetics that stabilize a water molecule within the β 3 subunit create a stable closed conformation of the binding site, inhibiting large-scale conformational changes [14].Assessing the structural rearrangements of the integrin binding site during RGD binding and evaluating the formation and disassembly of molecular contacts in a dynamic setting will provide insights for the future design of RGD-based drugs to inhibit conformational activation and prevent thrombosis.
Int. J. Mol.Sci.2024, 25, x FOR PEER REVIEW 2 of 17 increased thrombosis in some patients and exposure to neo-epitopes, causing an immune reaction that results in thrombocytopenia in others [13].RGD ligand mimetics that stabilize a water molecule within the β3 subunit create a stable closed conformation of the binding site, inhibiting large-scale conformational changes [14].Assessing the structural rearrangements of the integrin binding site during RGD binding and evaluating the formation and disassembly of molecular contacts in a dynamic setting will provide insights for the future design of RGD-based drugs to inhibit conformational activation and prevent thrombosis.The binding of an RGD peptide to αIIbβ3 integrin is a complex dynamic process involving structural changes of the ligand binding site (Figure 1E, F), accompanied by RGD motions and several sequential interactions of RGD with different integrin residues [17].RGD binding initiates molecular rearrangements of the αIIb β-propeller and β3 β-I domains forming the ligand binding site (Figure 1E, F).These rearrangements include a lateral movement of the backbone of β-Ser-123 in the β1-α1 loop, a downward shift of the α7 helix, a downward lateral motion of the carbonyl oxygen of β-Met-335 in the β6-α7 loop, and an increase in helicity in the α1 helix [17].These rearrangements change the position of the β-propeller domain, resulting in the opening of the headpiece and an increase in affinity for ligand binding [18][19][20][21].Unfortunately, experimental methods that can capture The binding of an RGD peptide to α IIb β 3 integrin is a complex dynamic process involving structural changes of the ligand binding site (Figure 1E,F), accompanied by RGD motions and several sequential interactions of RGD with different integrin residues [17].RGD binding initiates molecular rearrangements of the α IIb β-propeller and β 3 β-I domains forming the ligand binding site (Figure 1E,F).These rearrangements include a lateral movement of the backbone of β-Ser-123 in the β1-α1 loop, a downward shift of the α7 helix, a downward lateral motion of the carbonyl oxygen of β-Met-335 in the β6-α7 loop, and an increase in helicity in the α1 helix [17].These rearrangements change the position of the β-propeller domain, resulting in the opening of the headpiece and an increase in affinity for ligand binding [18][19][20][21].Unfortunately, experimental methods that can capture these rearrangements at high spatial and temporal resolutions are lacking, and computational modeling approaches present limitations in dynamic sampling.
In this study, we employed all-atom simulations using funnel metadynamics (FM) to enhance the sampling sufficiently to evaluate the dynamic interactions of an RGD peptide with the α IIb β-propeller and β 3 β-I domains of integrin.FM does not require a priori information about the interactions between RGD and integrin and incorporates an external history-dependent potential while applying a funnel-shaped restraint potential to limit RGD exploration of the unbound state.FM also provides structural flexibility for both protein and ligand and includes explicit water molecules for solvation and desolvation.It also considers conformational selection and induced-fit effects, which often play a role in identifying ligand binding modes [42].Our FM simulations showed that RGD binding to α IIb β 3 integrin triggers structural rearrangements within the β 3 β-I domain.Furthermore, our simulations provided a free-energy landscape with a low-energy binding mode surrounded by higher-energy prebinding states, demonstrating the existence of multiple binding modes.Our analysis also revealed several critical interaction residues distributed across the α IIb and β 3 subunits, aligning closely with findings from previous studies.In sum, this study demonstrates the effectiveness of funnel metadynamics in precisely evaluating the dynamics of RGD binding to α IIb β 3 integrin.

Structural Comparison between Cryo-EM and Reference Crystal Structures of Integrin α IIb β 3
The objective of this analysis was to evaluate the structural properties of the β-propeller and β-I domains of cryo-EM integrin α IIb β 3 used in this study by comparing these domains with the closed conformer (PDBID: 3T3P [15]) as a reference.The β-I domain of cryo-EM α IIb β 3 was structurally aligned to the β-I domain of the reference crystal, resulting in a Cα RMSD of 0.071 nm.The distances between specific backbone atoms in the cryo-EM α IIb β 3 and the corresponding atoms in the reference structure were measured: β-Met-335 oxygen (0.707 Å), α7 helix center of mass (COM) (0.955 Å), and β-Ser-123 nitrogen (0.915 Å).These comparisons confirm the cryo-EM structure's comparability to the reference crystal.

Structural Dynamics of Integrin β-I Domain: Insights from Funnel Metadynamics Simulations
Given that RGD binding to α IIb β 3 triggers molecular rearrangements in the β-I domain that propagate to the α subunit and to the lower legs for conformational activation [10,43,44], our objective was to evaluate the efficacy of funnel metadynamics (FM) simulations in capturing the mechanisms.Ten independent walkers were employed for FM, with those reaching 150 ns considered for structural analysis.Additionally, three 150 ns equilibrium molecular dynamics (MD), along with a 1 µs equilibrium MD simulation without RGD and divalent cations were run for comparison.
First, we examined the evolution of the helical content in the α1 helix.FM walkers with RGD and divalent cations displayed an increase in average helical content from ~70% to ~77% (Figure S1A,B).Equilibrium MD simulations on the same system showed variable results, with two out of three simulations exhibiting a decrease in helical content to ~53% and one showing an increase to ~76% (Figure S1C), suggesting that the presence of RGD and divalent cations can promote integrin activation but with lower probability than with biased ligand binding.Conversely, without divalent cations and RGD, the equilibrium simulation showed a decrease in average helical content to 63% (Figure S1D).
Next, we analyzed the dynamics of the α7 helix and the backbone oxygen atom of β-Met-335 in the β-I domain.We compared their distances from the corresponding positions in closed (PDB ID: 3T3P [15]) and open (PDB ID: 2VDR [16]) crystal structures.During FM simulations, these distances decreased, indicating a shift of the α7 helix and β-Met-335 from the bent towards the open state of integrin (Figures S2A,B and S3A,B).Equilibrium simulations on the same systems displayed more stable measures (Figures S2C and S3C), while the equilibrium simulation without divalent cations and RGD showed motion towards the closed crystal reference instead of the open (Figures S2D and S3D).
Overall, our analysis revealed that RGD binding and divalent cations are crucial for increasing helical rearrangements of the α1 helix and repositioning of the α7 helix and β-Met-335 away from the ligand binding interface of α IIb β 3 integrin in its bent closed conformation, towards their positions in the extended-open state.These findings affirm the capability of FM to capture molecular rearrangements of the β-I domain induced by RGD binding and the influence of divalent cations on these dynamics.

Assessing Sampling Efficiency: Evaluation of Funnel Metadynamics Simulations
To assess the efficacy of FM for simulations of the binding and unbinding between an RGD peptide and α IIb β 3 integrin, we conducted an evaluation of sampling efficiency across two key dimensions: the funnel space and the biased collective variable (CV) space (CV1 and CV2).Examination of the cumulative exploration of the ligand binding site for the RGD peptide, defined by the projections of the center of mass (COM) and contacts of RGD across all FM walkers, revealed a comprehensive exploration of both the funnel space and the biased CV space (Figure 2A,B).In the solvated state, where the RGD peptide exhibited no contact with the ligand binding site, the free energy surface appeared flat (refer to Figure 2C,D for CV1 = 3-5 nm).Conversely, in the bound state, a distinct energy well emerged, with the minimum corresponding to CV1 = 1.5 nm and CV2 = 0.1-0.3(Figure 2D).Through reweighted free energy analysis within the funnel space, we determined that the positions of the free energy minima were sufficiently distant from the funnel edges, suggesting that the funnel did not overly restrict the ligand within the binding site (Figure 2C).Overall, our findings demonstrate that FM offers robust sampling efficiency, ensuring a comprehensive exploration of the conformational landscapes across several free energy barriers.ure 2D).Through reweighted free energy analysis within the funnel space, we determined that the positions of the free energy minima were sufficiently distant from the funnel edges, suggesting that the funnel did not overly restrict the ligand within the binding site (Figure 2C).Overall, our findings demonstrate that FM offers robust sampling efficiency, ensuring a comprehensive exploration of the conformational landscapes across several free energy barriers.

Evaluation of RGD Interactions with the Binding Site of αIIbβ3 Integrin
Because the RGD peptide explored the integrin binding site within the funnel's cone, as well as the solvated state further from the cone, and did not experience bias from the funnel while in the bound state (Figure 2), we aimed to determine where binding occurred.Frames from FM simulations were extracted from all walkers where the COM of the RGD ligand was in close contact with the binding site around 1.5 nm on CV1, and close to 0 distance from reference contacts (CV2), and within 10 kJ/mol of the lowest energy point, (this region is indicated by a blue rectangle in Figure 2B,D).These extracted frames were used for all successive analyses.A contact between RGD and integrin was considered formed when a heavy atom within the RGD peptide and a heavy atom within a residue of the β-propeller or β-I domains were within 5 Å of each other [45].Many of the contacts reported by crystal structures and cryo-EM studies were confirmed in this
It has been previously reported that RGD contact with residues α-Asp-224 and α-Asp-232 is occasionally mediated by a water molecule, depending on the specific ligand and the conformation of integrin [16,17,22,25].To estimate the role of water in the integrin α-Asp-224 and α-Asp-232 contacts, the number of extracted contact frames where a single water molecule was detected within a specified distance (4, 3.3, and 2.6 Å) for both contacting residues was counted (Supplemental Table S1).For the α-Asp-224 and Arg-408 contact, the percentage of frames with water within 4, 3.3, and 2.6 Å were 55%, 40%, and 0.12%, respectively.For the α-Asp-232 and Arg-408 contact, the percentages were 100%, 99%, and 0.5%, respectively.Contacts reported to be exclusively mediated by water in crystal structure or cryo-EM studies but exhibiting direct contact within our criteria were β-Ala-218 (97%), β-Asp-251 (56%), α-Asp-163 (31%), and β-Tyr-166 (14%) (Figure 3C,D) [16,17,22,25].These contacts were also evaluated for water H-bond estimates (Supplemental Table S1).β-Ala-218 showed the most contact with Asp-410, with 84%, 43%, and 0% of the contact frames containing a water molecule within the respective distance criteria.β-Asp-251 showed 99%, 91%, and 5% of contact frames with a water molecule within the respective distances, and α-Asp-163 showed 100%, 98%, and 0.2% of contact frames with a water molecule within the respective distances.β-Tyr-166 showed the most contact with ligand residue Arg-408, with 99.5%, 88%, and 10% of contact frames showing a water molecule within the respective distances.β-Tyr-166 and β-Ala-218 showed contact with the other two ligand residues but to a much lesser extent (Supplemental Table S1).[15,16,27,29,30,[32][33][34].Red and blue squares indicate the labeled integrin residues, while the β-propeller and β-I domains are represented by a gray mesh and a pink mesh, respectively.The RGD peptide is shown in green, with atoms shown as small spheres (C,D).Contacts reported to be watermediated in previous structural studies [16,17,22,25] and confirmed here, with β-propeller and β-I domains represented in gray and pink, respectively, with a transparent surface.The RGD peptide is shown in green, with atoms shown as small spheres.The color scheme is as follows: oxygen in red, carbon in cyan, nitrogen in blue, hydrogen (on residues) or Ca 2+ (as a large sphere) in white, and Mg 2+ represented by an orange sphere.[15,16,27,29,30,[32][33][34].Red and blue squares indicate the labeled integrin residues, while the β-propeller and β-I domains are represented by a gray mesh and a pink mesh, respectively.The RGD peptide is shown in green, with atoms shown as small spheres (C,D).Contacts reported to be water-mediated in previous structural studies [16,17,22,25] and confirmed here, with β-propeller and β-I domains represented in gray and pink, respectively, with a transparent surface.The RGD peptide is shown in green, with atoms shown as small spheres.The color scheme is as follows: oxygen in red, carbon in cyan, nitrogen in blue, hydrogen (on residues) or Ca 2+ (as a large sphere) in white, and Mg 2+ represented by an orange sphere.

Table 1.
Comparison between α IIb β 3 integrin residues interacting with RGD in this study and previous investigations.The first column provides the list of residues from the two integrin subunits.The second column lists the percentage of frames capturing the contact.This percentage is calculated from the FM free energy surface as the total number of frames in which the RGD was found in the bound state.The third column provides references to previous studies that have demonstrated the existence of the contacts.* indicates residues with contacts for RGD that are water-mediated.

Discussion
The interplay between the α IIb and β 3 subunits in establishing the molecular contacts of integrin with an RGD peptide has been thoroughly explored using a variety of experimental and computational techniques, including X-ray crystallography, cryo-EM, MD simulations, crosslinking studies, loop swapping, and site-directed mutagenesis.While experimental methods provide invaluable insights, computational approaches offer higher spatial and temporal resolution, making them crucial for capturing detailed ligand-protein interactions, especially in a dynamic setting.Previous computational studies, however, have encountered limitations such as restricted sampling and the need for a priori information on integrin/RGD binding modes.These approaches have struggled to accurately capture the full transitions between unbound and bound states of the ligand, limiting their ability to identify complex binding modes.This study demonstrates the utility of funnel metadynamics in capturing structural changes at the binding site and exploring complex binding modes in a dynamic setting.
Application of FM to the integrin/RGD complex resulted in topological rearrangements within the β 3 subunit β-I domain (Figure S1), indicating that the transitions of the ligand between a solvated and a ligated state underlie integrin conformational activation.Consistent with this result, previous reports have shown that structural changes in the β-I domain correlate with the conformational transitions of integrin from inactive to active states [3,18,20].In particular, as ligand binding occurs and integrin moves from a bent-closed to an extendedopen conformation, the hybrid domain swings out, and within the β-I domain, the helicity of the α1 helix increases, the oxygen atom of β-Met-335 moves away from the binding site, and the α7 helix drops down [16,17,25,48], all supportive of our results.
Prior computational investigations relied on methodologies such as steered molecular dynamics (SMD), umbrella sampling (US), quantum mechanics (QM), and equilibrium molecular dynamics (MD).SMD simulations have been employed to extract the ligand from its binding site, facilitating the identification of contact-breaking events and the determination of the RGD unbinding force [28,29,31,32,36].Additionally, SMD has been utilized to drive an RGD peptide into the binding site, enabling the estimation of binding pathway energy, hydrogen bond formation, and the final distance between specific residues such as Asp-410 and the MIDAS Mg 2+ ion [36].Among the interactions investigated, those involving ligand residues Arg-408 and Asp-410, and α-Asp-224 and the MIDAS Mg 2+ ion, respectively, have been consistently identified as pivotal binding interactions within the α IIb β 3 integrin binding site.However, it is noteworthy that SMD simulations have typically employed pulling forces larger than physiological values.US simulations, on the other hand, have constrained the ligand's center of mass at various distances from the protein, allowing exploration of other degrees of freedom.Nonetheless, these simulations have not adequately distinguished the free energy differences attributable to these degrees of freedom.Equilibrium MD simulations have often initiated the ligand within the binding site; however, due to limited sampling, unbinding events have not been adequately captured [31,33,34].Conversely, starting the RGD ligand outside of the binding site has facilitated the observation of binding events but not unbinding [32].Quantum mechanical (QM) calculations have identified residues with the highest RGD interaction energies from the α IIb and β 3 chains, including residues such as α-Asp-224, α-Asp-159, β-Asn-215, and β-Lys-125, as well as the MIDAS Mg 2+ ion, and the Ca 2+ ions in the ADMIDAS and LIMBS sites [30].However, the inflexibility of the protein within the docking program, coupled with the computational complexity of QM calculations, has limited their ability to account for dynamic rearrangements of the target protein [38][39][40][41].None of these methods have fully captured the complete transitions of RGD between unbound and bound states, which are crucial for identifying complex binding modes.
To conclude, this study demonstrates the utility of funnel metadynamics in exploring structural rearrangements of the ligand binding site and integrin-RGD peptide interactions.It reveals molecular changes in the β 3 subunit of integrin upon RGD binding and provides a free-energy landscape with a low-energy binding mode surrounded by higher-energy prebinding states.Identifying both known and novel contacts enhances our understanding of ligand recognition mechanisms for α IIb β 3 integrin.The findings corroborate previous observations while providing additional molecular details, particularly regarding watermediated interactions and rarely observed contacts.Overall, the strong agreement between previous experimental and simulation data and our results highlights the reliability of FM as a method for studying complex protein/ligand systems.Funnel metadynamics could be applied in future research to inform therapeutic strategies targeting protein/ligand interaction processes.
To prepare the structure for simulation, we utilized the CHARMM-GUI Lehigh University, Bethlehem, PA, USA, Solution Builder [51,58].The setup included a pH of 7.0, employing the CHARMM36m force field [59], applying standard N and C-terminal group patching for all three protein chains, solvating with the CHARMM-modified TIP3P water model [60], and introducing 0.15 M Sodium Chloride (NaCl) to neutralize the charge of the system, using the Monte-Carlo placement method.The resulting system consisted of 162,986 atoms with 10,598 protein atoms, 50,696 water molecules, and 143 Cl, 150 Na, 6 Ca 2+ , and 1 Mg 2+ ions, arranged within a 120 × 120 × 120 Å box.Subsequently, the system underwent minimization and equilibration using the GRO-MACS 2022.5 software, Science for Life Laboratory, Stockholm, SE [61][62][63].Lennard-Jones interactions were truncated at 1.2 nm with the Verlet switching function, ranging from 1.0 to 1.2 nm.Long-range electrostatic interactions were calculated using the particle mesh Ewald (PME) summation method with a 0.12 nm grid spacing [64].The systems were minimized using the steepest descent algorithm for 5000 steps or until a maximum force of 1000 kJ mol −1 nm −1 was reached, with a step size of 0.01 nm.Positional restraints for both minimization and equilibration were set at 400 kJ mol −1 nm −2 for protein backbone atoms and 40 kJ mol −1 nm −2 for protein side chain atoms.Subsequently, the temperature was raised to 310 K using velocity rescaling with the Nose-Hoover thermostat and equilibrated for 125,000 steps with 0.001 ps timesteps (125 ps total), implementing the canonical ensemble (NVT), which maintains a constant number of atoms, volume, and temperature.Subsequently, positional restraints were released, and production runs were conducted in the isothermal-isobaric ensemble (NPT), which maintains a constant pressure of 1 atm using isotropic Parrinello-Rahman coupling [65] and a temperature of 310 K using the Nose-Hoover thermostat.Covalent bonds involving hydrogen atoms were constrained with the LINCS algorithm [66], allowing for a timestep of 0.002 ps.This equilibrated system was used to produce three replicas of equilibrium simulations of 150 ns each.Additionally, to facilitate steered molecular dynamics (SMD) and subsequent funnel metadynamics (FM) simulations, a restraint based on a harmonic potential with a spring constant of 10 kJ/mol was applied to the alpha carbon (Cα) atoms in residues 96 and 107 within the β-propeller domain.These residues were chosen based on the evaluation of root mean squared fluctuation (RMSF) of individual integrin Cα atoms over the first 50 ns of equilibrium MD simulations, which indicated that these atoms exhibit the lowest fluctuation across the three replicas.

Steered Molecular Dynamics Simulations
The MOVINGRESTRAINT bias, implemented in PLUMED 2.8.2, was utilized to extract the RGD ligand from the binding site while maintaining the funnel restraint.PLUMED is an open-source, community-developed library [67].The resulting trajectory served as the basis for initializing the walkers.The restraint was initiated near the ligand's initial center of mass (COM) position at 0.4 nm and gradually extended to 4.5 nm over 2,500,000 steps (5 ns) with a linear increase in the spring constant during this period.Initially set at 0 kJ/mol, the spring constant was progressively raised to 1000 kJ/mol at the 5 ns mark before reverting to 0 kJ/mol at the 10 ns mark.The total simulation time was 7 ns.Frames and PLUMED output were saved every 500 steps (1 ps).Further details on the funnel setup are elaborated in the subsequent section.

Funnel Metadynamics Simulations
Following the Steered Molecular Dynamics (SMD) protocol, a series of multiple walker well-tempered Funnel Metadynamics (FM) simulations were conducted using GROMACS 2022.5, integrated with PLUMED 2 [68].This approach employs multiple walkers with a "cone and cylinder" or funnel-shaped restraint, in conjunction with well-tempered metadynamics, aimed at enhancing the sampling within the target ligand binding site and its solvated state.The funnel constraint was imposed solely when the ligand approached the edge of the funnel, as depicted in Figure 4A,B.Simultaneously, the cylinder allowed controlled exploration of the solvated state, thereby increasing the sampling of RGD in both bound and solvated states.To prevent the RGD ligand from diffusing into the bulk solvent, a restraint in the form of a harmonic upper wall was positioned at the end of the cylinder, Two collective variables (CVs) were selected to drive the metadynamics bias: CV1, representing the center of mass (COM) position of the RGD peptide along the funnel's zaxis, and CV2, measuring the distance from a reference contact map (unitless) (Figure 4B,  C).These contacts involved residues from the αIIb subunit that interact with Arg-408 (α-Tyr-189, α-Tyr-190, α-Asp-224), residues from the β3 subunit that interact with Asp-410 (β-Asn-215, β-Tyr-122, β-Ser-123), and the Mg 2+ ion that interacts with Asp-410 (Figure 4C).
CV2 called the CONTACTMAP within PLUMED, was computed using a rational switching function, defined as: The anchoring point and origin of the funnel were defined as atom number 7655, corresponding to a Cα atom on β-Pro-163 in the β-I domain, thereby centering the cone on the α IIb β 3 binding site.The second funnel point was chosen on the docked ligand to ensure alignment of the α IIb β 3 binding site and bound ligand within the funnel.Specific funnel parameters were set as follows: a cylinder radius of 1.0 Å, an amplitude of 0.35 radians for the cone section, and a switching point between the cone and the cylinder set at 2.0 nm.The center of mass (COM) calculation of the ligand's position along the z-axis of the funnel utilized all heavy atoms within the RGD ligand chain (residue 408 to 410).
Multiple walker metadynamics allowed concurrent contributions from multiple simulations (walkers) to the metadynamics bias, thereby enhancing sampling efficiency within a shorter real-time duration.Ten walkers were initiated from points within the SMD trajectory, each commencing at a frame nearest to the beginning of the SMD simulation and spaced at regular 0.3 nm intervals along the funnel, as illustrated in Figure 4B, ranging from 1.6 nm to 4.1 nm.Walker 5, initiated at 3.1 nm, attained a simulation time of 121 ns, revealing an error in calculating the funnel grid, potentially stemming from inadequate restraints confining the ligand within the funnel boundaries or adjustments in the funnel leading to ligand placement outside the funnel grid between timesteps.Each of these walkers, except where noted above, reached 150 ns.
CV2 called the CONTACTMAP within PLUMED, was computed using a rational switching function, defined as: where d 0 = 0.0, n = 6, m = 2n, r 0 = 1 nm is the switch distance, and r is the distance of the atom from its contact pair.Each switched distance is then squared and summed, and the resulting value is then subtracted from the same calculation with the switched reference values.This equation, known as the CMDIST function within PLUMED, is written as: where S ref1,ref2,ref3,. . .are the switched output of the reference contacts and S 1,2,3. . .are the switched output of the same contacts from the current timestep.The first CV allows the bias to further enhance sampling of the solvated and bound state, and the second CV was chosen to distinguish different binding modes of the ligand in the binding site, if present.The Gaussian deposition rate was every 500 steps (1 ps) with a height of 2.0 kJ/mol and a width of 0.1 nm for CV1 and 0.5 for CV2, based on equilibrium simulations.The bias factor was set to 20.Additionally, the WHOLEMOLECULES feature was employed in PLUMED 2 for all protein atoms to facilitate the reconstruction of these molecules within PLUMED in case of crossing the periodic boundary, thereby preventing calculation artifacts.Output from PLUMED and coordinates (frames) from GROMACS were saved every 500 steps (1 ps).

Equilibrium Molecular Dynamics Simulations without Divalent Cations and RGD
The full-length α IIb β 3 integrin structure from [2], modified to match the canonical sequence at uniprot.org, was embedded into a membrane and solvated at pH 7.0 using CHARMM-GUI Membrane Builder [51].To align the protein with the membrane plane, the Orientation Option "Run PPM 2.0" was employed, taking input from both integrin subunits.A lipid ratio of 7:3 POPE:POPC (palmitoyl-oleoyl-phosphatidylethanolamine:palmitoyloleoyl-phosphatidylcholine) was utilized, resulting in a membrane configuration similar to cellular membranes.Sodium chloride at a concentration of 0.15 M was incorporated via the distance ion placing method.The resulting system comprised 359,339 atoms, with 540 lipid molecules, 87,641 water molecules, 292 Na ions, and 233 Cl ions, arranged within a box of dimensions 132 × 132 × 217 Å.Note that no divalent cations or RGD ligands were included.
The system underwent minimization using the steepest descent algorithm for 5000 steps or until a maximum force of 1000 kJ mol −1 nm −1 was reached, with a step

Figure 1 .
Figure 1.All-atom representation of αIIbβ3 integrin.(A,B) Different viewpoints of full-length integrin in quick surf representation.The αIIb subunit is in gray; the β3 subunit is in pink.(C) The ligand binding site at the interface between αIIb β-propeller and β3 β-I domains (quick surf representation) with a bound RGD peptide (licorice representation).(D) The ligand binding site with the representation in C is made transparent and overlaid with integrin (new cartoon representation).(E) Overlay of the closed cryo-EM structure (Cyan) used in this study with the reference closed (Orange, PDBID: 3t3p [15]) and open (Green, PDBID: 2VDR [16]) structures.Specific atoms and regions are highlighted within the β-I domain of the β3 chain.These structures are known to show the biggest changes between closed and open crystal structures.(F) zoomed in view of the β-I domain in E. Gray and pink represent the αIIb and β3 subunits, respectively.In the ligand (C,D), cyan = carbon, blue = nitrogen, red-oxygen.(D,E) white spheres = Ca 2 + , orange sphere = Mg 2+ .

Figure 1 .
Figure 1.All-atom representation of α IIb β 3 integrin.(A,B) Different viewpoints of full-length integrin in quick surf representation.The α IIb subunit is in gray; the β 3 subunit is in pink.(C) The ligand binding site at the interface between α IIb β-propeller and β 3 β-I domains (quick surf representation) with a bound RGD peptide (licorice representation).(D) The ligand binding site with the representation in C is made transparent and overlaid with integrin (new cartoon representation).(E) Overlay of the closed cryo-EM structure (Cyan) used in this study with the reference closed (Orange, PDBID: 3t3p [15]) and open (Green, PDBID: 2VDR [16]) structures.Specific atoms and regions are highlighted within the β-I domain of the β 3 chain.These structures are known to show the biggest changes between closed and open crystal structures.(F) zoomed in view of the β-I domain in E. Gray and pink represent the α IIb and β 3 subunits, respectively.In the ligand (C,D), cyan = carbon, blue = nitrogen, red-oxygen.(D,E) white spheres = Ca 2+ , orange sphere = Mg 2+ .

Figure 2 .
Figure 2. Exploration of the RGD peptide of the funnel and biased CV space during FM simulations.(A) RGD center of mass (COM) exploration within the funnel space.(B) RGD exploration in the space of the biased collective variables (CVs), with CV1 representing RGD COM along the funnel zaxis and CV2 indicating the difference in contact maps.(C) Two-dimensional projection of the reweighted free energy surface of RGD COM within the funnel space.(D) Three-dimensional representation of RGD energy as a function of the biased collective variables in (C).The blue squares in panels (B,D) denote the minimum of the free energy surface, which is used to determine CV space boundaries for frame extraction and contact analysis.The arrow indicates the direction of decreasing energy: lighter colors correspond to higher energy values; darker colors correspond to lower energy values.

Figure 2 .
Figure 2. Exploration of the RGD peptide of the funnel and biased CV space during FM simulations.(A) RGD center of mass (COM) exploration within the funnel space.(B) RGD exploration in the space of the biased collective variables (CVs), with CV1 representing RGD COM along the funnel z-axis and CV2 indicating the difference in contact maps.(C) Two-dimensional projection of the reweighted free energy surface of RGD COM within the funnel space.(D) Three-dimensional representation of RGD energy as a function of the biased collective variables in (C).The blue squares in panels (B,D) denote the minimum of the free energy surface, which is used to determine CV space boundaries for frame extraction and contact analysis.The arrow indicates the direction of decreasing energy: lighter colors correspond to higher energy values; darker colors correspond to lower energy values.

Figure 3 .
Figure 3. Structural analysis of the interactions between the RGD peptide and the ligand binding site of αIIbβ3 integrin.(A) Residues in the αIIb domain and (B) in the β-I domain that interact with the RGD peptide, as identified in this study and consistent with crystal structure and cryo-EM studies[15,16,27,29,30,[32][33][34].Red and blue squares indicate the labeled integrin residues, while the β-propeller and β-I domains are represented by a gray mesh and a pink mesh, respectively.The RGD peptide is shown in green, with atoms shown as small spheres (C,D).Contacts reported to be watermediated in previous structural studies[16,17,22,25] and confirmed here, with β-propeller and β-I domains represented in gray and pink, respectively, with a transparent surface.The RGD peptide is shown in green, with atoms shown as small spheres.The color scheme is as follows: oxygen in red, carbon in cyan, nitrogen in blue, hydrogen (on residues) or Ca 2+ (as a large sphere) in white, and Mg 2+ represented by an orange sphere.

Figure 3 .
Figure 3. Structural analysis of the interactions between the RGD peptide and the ligand binding site of α IIb β 3 integrin.(A) Residues in the α IIb domain and (B) in the β-I domain that interact with the RGD peptide, as identified in this study and consistent with crystal structure and cryo-EM studies[15,16,27,29,30,[32][33][34].Red and blue squares indicate the labeled integrin residues, while the β-propeller and β-I domains are represented by a gray mesh and a pink mesh, respectively.The RGD peptide is shown in green, with atoms shown as small spheres (C,D).Contacts reported to be water-mediated in previous structural studies[16,17,22,25] and confirmed here, with β-propeller and β-I domains represented in gray and pink, respectively, with a transparent surface.The RGD peptide is shown in green, with atoms shown as small spheres.The color scheme is as follows: oxygen in red, carbon in cyan, nitrogen in blue, hydrogen (on residues) or Ca 2+ (as a large sphere) in white, and Mg 2+ represented by an orange sphere.

Figure 4 .
Figure 4. Representation of the integrin/ligand complex including the funnel structure, and indication of the collective variables for FM simulations.(A) New cartoon representation of the αIIb βpropeller domain (gray) and the β3 β-I domain (pink), forming the ligand binding site.The RGD peptide is in licorice representation, while the funnel is shown in transparent orange.The first collective variable (CV1), corresponding to the funnel z-axis, is shown with a red line.(B) Initial spatial positions of the RGD peptide for the various simulation walkers along CV1.The black sphere denotes the center of mass (COM) of the RGD peptide along the funnel z-axis, as derived from SMD, serving as the starting coordinates for the RGD ligand for each of the ten walkers.(C) Contacts considered in the calculation of the second collective variable (CV2).These contacts are formed between residues from the αIIb subunit that interact with Arg-408 (α-Tyr-189, α-Tyr-190, α-Asp-224), residues from the β3 subunit that interact with Asp-410 (β-Asn-215, β-Tyr-122, β-Ser-123), and the Mg 2+ ion that interacts with Asp-410.Dashed lines represent the eight contacts utilized to measure the deviation, in distance, between the initial reference contacts and those observed throughout the simulation.Note that β-Asn-215 forms two interactions with Aps-410.Atoms are color-coded as follows: oxygen (red), carbon (cyan), nitrogen (blue), calcium (white spheres), magnesium (orange sphere), and center of mass of atoms forming contacts with RGD (black spheres).

Figure 4 .
Figure 4. Representation of the integrin/ligand complex including the funnel structure, and indication of the collective variables for FM simulations.(A) New cartoon representation of the α IIb β-propeller domain (gray) and the β 3 β-I domain (pink), forming the ligand binding site.The RGD peptide is in licorice representation, while the funnel is shown in transparent orange.The first collective variable (CV1), corresponding to the funnel z-axis, is shown with a red line.(B) Initial spatial positions of the RGD peptide for the various simulation walkers along CV1.The black sphere denotes the center of mass (COM) of the RGD peptide along the funnel z-axis, as derived from SMD, serving as the starting coordinates for the RGD ligand for each of the ten walkers.(C) Contacts considered in the calculation of the second collective variable (CV2).These contacts are formed between residues from the α IIb subunit that interact with Arg-408 (α-Tyr-189, α-Tyr-190, α-Asp-224), residues from the β 3 subunit that interact with Asp-410 (β-Asn-215, β-Tyr-122, β-Ser-123), and the Mg 2+ ion that interacts with Asp-410.Dashed lines represent the eight contacts utilized to measure the deviation, in distance, between the initial reference contacts and those observed throughout the simulation.Note that β-Asn-215 forms two interactions with Aps-410.Atoms are color-coded as follows: oxygen (red), carbon (cyan), nitrogen (blue), calcium (white spheres), magnesium (orange sphere), and center of mass of atoms forming contacts with RGD (black spheres).