Single molecule force spectroscopy data and BD- and MD simulations on the blood protein von Willebrand factor

We here give information for a deeper understanding of single molecule force spectroscopy (SMFS) data through the example of the blood protein von Willebrand factor (VWF). It is also shown, how fitting of rupture forces versus loading rate profiles in the molecular dynamics (MD) loading-rate range can be used to demonstrate the qualitative agreement between SMFS and MD simulations. The recently developed model by Bullerjahn, Sturm, and Kroy (BSK) was used for this demonstration. Further, Brownian dynamics (BD) simulations, which can be utilized to estimate the lifetimes of intramolecular VWF interactions under physiological shear, are described. For interpretation and discussion of the methods and data presented here, we would like to directly point the reader to the related research paper, “Mutual A domain interactions in the force sensing protein von Willebrand Factor” (Posch et al., 2016) [1].

methods and data presented here, we would like to directly point the reader to the related research paper, "Mutual A domain interactions in the force sensing protein von Willebrand Factor" (Posch et al., 2016) [1].
& This data in brief article describes how different methods, such as AFM based single molecule force spectroscopy and BD-and MD simulations can be used to study the protein VWF. Single molecule force spectroscopy (SMFS), which is described in this article, is a technique to nondestructively probe the forces and the dynamics between molecules under physiological conditions. In combination with our chemical expertize, this powerful tool is applicable for almost every receptor/ligand combination, which offers many possibilities for new collaborations.
The detailed description of the calculation of tensile force profiles along a protein chain from BD simulations might be a guideline for other researcher pursuing similar interests.
The fitting of rupture forces versus loading rate profiles in the MD loading-rate range show an elegant way for the comparison of data from experimental AFM experiments and theoretical MD simulations and provide for collaboration between experimentalists and theorists.

Data
In Fig. 1 we provide deeper understanding of the recording of force distance cycles and advanced data evaluation. The process of specifity proof measurements using the example of VWF A1/A2 interaction studies is discussed and a typical loading rate dependence plot is shown in Figs. 2 and 3.
The BD simulation part ((Eqs. 1-6) and Fig. 4) gives a detailed description of the calculation of force profiles for protein chains and the fitting of rupture forces versus loading rate profiles in the MD loading-rate range can be found in Fig. 5   showing a specific VWF A1/A2 unbinding event, which was marked with a polynomial fit. The unbinding force was determined from the jump at the point of dissociation as indicated. In the blocking experiment ( Fig. 2) no specific unbinding events occurred (Insert in A). At least 1000 FDCs were recorded for each loading rate to create a typical force distribution as shown in (B). The distribution was taken from the A1/A2 interaction measurement for a velocity of 400 nm/s. A Gaussian distribution was fitted to the first peak and forces within the interval m 7 σ were used for further analysis in the loading rate dependence plot (Fig. 3). The loading rate r was calculated by multiplying the pulling velocity v with the effective spring k eff (slope at the rupture) of the system.   utilized. The actual spring constant was determined using the thermal noise method [2]. AFM cantilever tips as well as mica supports were chemically functionalized with different VWF constructs [1].

Materials
All chemicals were used in the highest available purity. 3-Aminopropyl-triethoxy silane (APTES; Sigma Aldrich, Vienna, Austria) was distilled at low pressure and stored under argon in sealed crimp vials over dry silica gel (to avoid polymerization) at À20°C. MilliQ (Millipore, Massachusetts, USA) purified water was used for all aqueous solutions. Triethylamine (TEA, Sigma Aldrich, Vienna, Austria) was stored under argon and in the dark to avoid amine oxidation. Chloroform was purchased from J.T. Baker (Griesheim, Germany), Argon and N 2 from Linde Gas GmbH (Stadl-Paura, Austria). Concentrated HCl was purchased from Sigma Aldrich (Vienna, Austria). Ethylenediaminetetraacetic acid (EDTA) and Tris base were purchased from VWR International (Vienna, Austria), Hepes and NiCl 2 were obtained from Merck (Darmstadt Germany) and TCEP (tris(2-carboxyethyl)phosphine) hydrochloride from Invitrogen (Vienna, Austria). Disulfide-tris-NTA was generously provided by the Tampé lab, Biocenter, Frankfurt am Main, Germany. The heterobifunctional crosslinker maleimide-PEG 27 -NHS was purchased from Polypure (Oslo, Norway). N-succinimidyl 4-(dimethoxymethyl)benzoate (SDMB) was synthesized as described in [3]. The cDNAs' coding for recombinant human VWF constructs containing the A1A2 (aa 1230-1672) construct, the single A1 (aa 1230-1463) and the single A2 (aa 1494-1672) were cloned into the mammalian expression vector pIRES neo2 [4]. All VWF constructs were labeled with a His 6 -tag. Mutations were inserted by site-directed mutagenesis employing the QuickChange kit (Stratagene). The vectors were used to transform Top10 super competent cells (Invitrogen) and sequenced. Plasmid purification was performed using the Endofree Plasmid Maxi Kit (QIAGEN). HEK293 cells were cultured in Dulbecco Modified Eagle Medium (DMEM, Invitrogen) with 10% [v/v] fetal bovine serum (Invitrogen) and 1% penicillin/streptavidin at 37°C and 5% CO 2 . These cells were transfected with the VWF vectors using Lipofectamine 2000 (Invitrogen) according to the manufacturer's instructions and recombinant expression of VWF variants was performed as previously described [5]. His-tagged VWF domain constructs were purified employing the His-Pur Ni-NTA Resin (Thermo Fisher Scientific) according to the manufacturer's instruction for purification of His-tagged proteins using a gravity-flow column. Mica sheets were bought from Christine Groepl, Electron Microscopy (Tulln, Austria). For aqueous solutions, TBS buffer (50 mM Tris, 150 mM NaCl, pH 7.4) and Hepes buffer (prepared from a 1 M solution of Hepes acid by adjustment of pH 7.5 or pH 9.6as stated in the textwith 20% NaOH) were used.

Methods
Interactions were probed by conducting force-distance-cycles (FDCs) at different loading rates r, i.e., at nine different velocities (50, 100, 200, 400, 600, 800, 1200, 2000 and 3000 nm/s). To gain reliable statistics, at least 1000 FDC were recorded at each pulling speed. The position of the tip relative to the surface was changed every 200 FDC, so as to statistically avoid position dependent artifacts. Unbinding events in the recorded FDC were marked with a polynomial fit (Fig. 1(A)). Binding events could be discerned from nonspecific adhesion by a characteristic parabolic force signal due to the elastic properties of the linkers. To additionally prove the specificity of the interactions, blocking experiments were performed. In this process, the ligand on the tip (e.g. VWF domain A1) was incubated with free receptors (e.g. VWF domains A2) in solution for two hours (Fig. 2). The free receptors saturate the ligand on the tip, and thus block the interaction between the ligand and the receptor on the sample surface. Thereafter, almost no specific interactions were observed in the FDC and the binding probability (BP), defined as the ratio between FDCs showing an unbinding event relative to the total number of FDCs, dramatically decreased (Fig. 2).
Probability density functions (pdf) of the measured forces were computed for each loading rate. A Gaussian distribution was fitted to the first peak of the pdf and forces within the interval m 7σ were used for further analysis (Fig. 1(B)). The unbinding forces were plotted against the logarithm of the loading rate r, which is given by the product of the effective spring constant k eff (slope at the rupture of an unbinding event, see Fig. 1(A)) and the pulling velocity v. The data points in this loading rate dependence plot (LRD) represent a single unbinding event each and were fitted with a single energy barrier binding model [6], using the maximum likelihood fitting routine [7] to obtain the kinetic off-rate k off and the width of the energy barrier x β . A typical example for such a data cloud and its fit is shown in Fig. 3.

Calculating tensile force profiles (under shear flow conditions) from BD simulations
We used Brownian dynamics including hydrodynamic interactions [8] to simulate the behavior of VWF under shear flow conditions. The mean tensile force along the backbone at the grafted chain end was then plugged into the Bell Evans Eq. (1), see also [1]: to gain an estimate for the lifetime of the A1/A2 complex under physiological conditions. In the following, we will expand on the model and the parameter values we used and how we calculated the tensile force profiles from the simulation.

The model
A linear chain of VWF protomers is modeled as N beads interacting by a Lennard-Jones potential. The backbone of the chain is additionally held together by a spring interaction. The potential energy of the chain therefore reads where r i is the position vector of bead i and r ij is defined by the norm of the distance vector r i À r j . The parameters ε; a and κ are the cohesion strength, the bead radius and the spring constant, respectively. The repulsive potential of the surface was chosen to be: in accordance to [13].

Integration scheme
A detailed description of the integration scheme can be found in ( [14], Eq. (18)) and ( [13], Eq. (1)). A simplified notation of the algorithm is given by this formula where v shear is the velocity of bead i due to the shear flow, z i is the third component of r i , μ 0 is a mobility given by Stoke's formula, μ ij is the 3 dimensional submatrix of the Rotne-Prager-Blake tensor μ with μ 3i À 2;3j À 2 being the upper left entry of that submatrix as defined in Eq. (5), F ij is the force, that bead j exerts on the bead i andx andẑ are dimensionless unit vectors in x and z direction, respectively. The terms on the right hand side of Eq. (4) account for (in order left to right): the shear flow, all deterministic forces acting on bead i, a correction due to the spatial dependence of the random velocity, the random velocity (assumed to fulfill the fluctuation-dissipation theorem).
Hydrodynamic interactions are treated through the mobility matrix and a detailed description of the employed mobility matrix can be found in [14] Eq. (7)ff.

Parameter values
In our coarse-grained bead spring model (Eq. 2) each bead corresponds to a protomer of VWF with a reported radius of gyration of a¼ 30 nm [9]. The cohesion strength was chosen to be 2k B T, since this value leads to a collapsed conformation of the VWF multimer in the untethered case [10] and thus resembles the physiological situation of VWF in the bloodstream The spring constant for the interaction between two adjacent beads was set to κ ¼ 20 3 k B T nm . The first bead was put at the no slip boundary and fixed throughout the course of the simulation. We added a weak repulsive potential with a short range of σ R ¼45 nm to make sure beads could not penetrate the wall. All other bead's positions were updated using a time step of Δt ¼ 55 ns over the course of 10 8 steps covering a time frame of 5.5 s. The time step Δt was calculated assuming T to be room temperature and the dynamic viscosity of the surrounding fluid to be 0.89 mPa s. Prior to the simulation run, an equilibration run of at least 10 5 steps has been performed.

Tensile force calculation
Block averages of the distance of adjacent beads r i;i þ 1 were recorded over 100 steps. Those were again averaged at the end of the simulation run, yielding a mean distance 〈r i;i þ 1 〉 of adjacent beads.
Mean tensile forces f i where then calculated from the equation 2.3. Fitting of rupture forces versus loading rate profiles in the MD loading-rate range Rupture forces were fitted in a broad loading-rate range covering both AFM and MD regimes, using the model recently developed by Bullerjahn, Sturm, and Kroy (BSK) [11] (Fig. 4(C) and (D), [1]). In this model, three parameters were used for the fitting: the energy barrier to overcome during rupture (E), the distance from the bound to the transition rupture state (x b ), and the diffusion constant of the system along the pulling reaction coordinate D. For a given loading rate LR and parameters (E, xb, D), the BSK method yields the probability distribution of rupture forces: P(F; LR, E, xb, D). Accordingly, a maximum likelihood approach was used in conjunction to the BSK model to determine the set of parameters that best described the entire set of rupture forces recovered at different applied loading rates. Fitting was carried out considering, the AFM data of the VWF A1-A2 complex either in its wildtype (A1/A2) or its bridged (A1/[A2]) form (Fig. 5). AFM data used for fitting was obtained as explained as follows. Histograms of the AFM rupture forces were computed for each pulling velocity. Peaks in the histograms were extracted and a Gaussian distribution was fitted around the peak with the lowest force. This peak was considered for the analysis because it corresponds in most of the cases to the fingerprint of the VWF A1-A2 wild-type rupture process (see Fig. 1(B)). Resulting Gaussian distributions were considered as the input for the maximum likelihood BSK-based approach. The loading rate (LR) was estimated as the product LR ¼V oke 4, where V is the pulling velocity and ke is the effective spring constant of the linkers, cantilever, and VWF proteins. ke was determined from the slope of each force distance cycle. o 4 denotes average over cycles. Note that the elastic constant ks of the pulling apparatus (cantilever and linkers) should be used instead. Nevertheless, if the protein is sufficiently stiff and thereby its elastic constant is substantially large (kp-1), then the effective elastic constant corresponds to that of the pulling apparatus (1/ke ¼1/ks þ 1/kpE 1/ks). Indeed, kp was found to almost three orders of magnitude larger than ke, thus justifying our use of ke as an estimate of ks. Fitting was also performed by combining AFM data with MD data, for both A1-A2 wildtype and bridged construct variants. MD rupture forces recovered from the smoothed force-distance profiles were used (Fig. 4, [1]). For the wild-type construct, all MD rupture forces were considered. For the bridged construct, the largest force value was excluded as it was obtained from a simulation involving the unfolding of A2, a situation that is prevented for this construct. The pulling velocity in the force-probe MD simulations was V¼0.2 m/s and the spring constant of the two virtual springs (connected in series) was 415.145 pN/nm. Therefore the MD loading-rate was found to be LR ¼8.3 Â 1010 pN/s. The activation rate k0 in the absence of force was estimated as a function of the three fitting parameters based on Kramers rate theory as [12]: where k B is the Boltzmann temperature and T the temperature.