A mechano-reactive coarse-grained model of the blood-clotting agent von Willebrand factor

When humans bleed, a large blood glycoprotein, von Willebrand Factor (vWF), plays a vital role in the hemostasis process in a rapid flow field.^1–3 Specifically, vWF contributes to the function of aggregating platelets to collagen around a damaged vessel wall. The plasma vWF monomer contains the sequential domains of D′D3-A1-A2-A3-D4-C1-C2-C3-C4-C5-C6-CT/CK⁴ [for available structural information, see PDB 1M10 (A1), PDB 3GXB (A2), and PDB 4DMU (A3) in Protein Data Bank^5–7]. Both A1 and A3 domains are claimed to be involved in the hemostasis process through their binding to platelets and collagen under rapid flow conditions (i.e., for shear rates above 1500 s⁻¹).^8,9 Nonetheless, it has recently been noted that the binding between A1 and GPIbα (of the GPIb-IX-V complex, a profuse membrane receptor complex) to mediate platelet adhesion could be prohibited due to the adjacent, folded A2 domain.^10,11 An increased hydrodynamic field may induce unfolding of the A2 domain to allow A1–platelet and A3–collagen interactions.¹² The native multimeric glycoprotein, ultralarge vWF (ULVWF) is sliced into relatively small fractions by the plasma protease ADAMTS13 (ADAMTS13: A Disintegrin and Metalloprotease with a ThromboSpondin type 1 motif, member 13)¹³ to prevent uncontrolled aggregation of platelets. This size-regulation process can only be functional when the proteolytic site buried in A2 is exposed upon A2 domain unfolding.^14–16 There are at least two kinds of diseases—type 2 von Willebrand Disease (vWD) and thrombotic thrombocytopenic purpura (TTP)—related to dysfunctional vWF. Some vWD cases may be related to insufficient resistance to A2 domain unfolding, resulting in excessive molecular scission and thus insufficient vWF length to stop bleeding.^17,18 On the contrary, TTP is caused by excessive circulating uncleaved ULVWF, i.e., failure in the size regulation process that may have dysfunctional A2 domain behavior as part of its cause.^13,19

A distinction of vWF polymer transitions at two separate length scales (i.e., multimer vs domain scale) was previously observed. Particularly, there were two different conformational changes important to vWF functionality;²⁰ one is the whole polymer unraveling, and the other is the individual A2 domain unfolding. For the former process, Schneider et al. found that vWF polymers adopted a globular conformation in quiescent or low shear flow conditions, then transitioned to an unraveled conformation when the shear rate was sufficiently high (usually over 1000 s⁻¹).²¹ Previous simulations of vWF multimers adopted coarse-grained (CG) models with a single spherical bead representing the entire monomer.^22–24 In those works, each bead represented a protomer, and adjacent protomers were connected by relatively rigid harmonic springs.²¹ Such models demonstrated the conformational changes in whole vWF polymer under various shear rate conditions but could not directly simulate inner monomer unfolding. Breakthrough experimental analyses of inner monomer unfolding were made by Zhang et al.²⁵ Using single molecule force spectroscopy, they characterized the mechanical response of the inner structure of a vWF monomer. According to this study, the force-sensitive domain is A2, which resides in the central region of each monomer. Also, the A2 domain showed a nontrivial extension when unfolding in response to the external force.²⁵ Moreover, they elucidated that A2 domain unfolding events are characterized by a significant force decrease at constant extension.²² Based on the evidence advanced from Zhang et al., many researchers applied molecular dynamics (MD) studies to seek a deeper understanding of A2 unfolding events at the atomic level. Chen et al. and Dong et al. expressed the A2 domain unfolding to be a multistate process through all-atom MD simulations.^26,27 Like the process suggested in experiments, these all-atom studies applied external forces to the termini of isolated A2 domain to move them apart with a constant pulling speed. Multistate unfolding events of A2 domain were observed with a force drop occurring for each intermediate unfolding event. Although these studies are able to show detailed insights of the A2 unfolding phenomenon, the lower limit of pulling speed that can be used in MD simulations due to the simulation time is still three orders of magnitude higher than the upper capability in experiments, which could lead to unrealistic high unfolding forces.²⁷ In addition, because the vWF monomer crystal structure has not yet been fully revealed, all-atom MD simulations are not possible for exploring multimer conformation under flow conditions. Previously, we also studied shear induced A2 domain force excursions by using single finitely extensible nonlinear elastic (FENE) springs to represent A2 domains in a vWF multimer CG model.^12,28,29 This CG model has the ability to elucidate the conformational change in vWF multimers and A2 force excursion with a relatively low computational cost. However, the A2 domain unfolding phenomenon, including the dramatic force decrease that occurs upon unfolding, could still not be observed. To capture vWF multimer conformation and A2 domain unfolding behavior at the same time, a new CG model is needed. It is also important to point out that the process of ADAMTS13 cleavage is much more complicated than A2 domain unfolding; nonetheless, it is still meaningful to determine the unfolding behavior of A2 because it is believed to be necessary for cleavage. Despite this, there is no CG model reported in the literature that explicitly describes the A2 domain unfolding process.

In this work, we introduce a new activated A2 monomer model (AA2MM) to investigate A2 domain unfolding. AA2MM is a CG model that employs spring systems to mimic the mechanical force-extension reaction under tensile loading (e.g., due to hydrodynamic forces in flow). AA2MM was parameterized at the monomer level using the data from Zhang et al.²⁵ and then also validated at the polymer level using the data from Fu et al.³⁰ Finally, vWF polymer unfolding and A2 domain unfolding at the monomer scale were examined (or predicted) for bulk shear and elongational flow.

In order to investigate the mechanical behavior of the A2 domain unfolding process under shear flow conditions, Ouyang et al. introduced a CG model of the vWF monomer that consisted of two beads connected by a finitely extensible nonlinear elastic (FENE) spring; the FENE spring was used to represent the force-extension property of A2 between the A1 and A3 domains in vWF.²⁹ In our prior model and the model presented here, each bead was assumed to represent all the domains on either side of the A2 domain. Thus, we refer to the bead containing the A1 domain as A1 and similarly, the bead containing the A3 domain as A3. The two spherical beads are of the same radius^3,20,31 and are connected by a FENE spring (i.e., the model A2 domain) to represent a monomer. To form polymers, relatively stiff harmonic springs were employed to connect adjacent monomers in a head-to-head, tail-to-tail manner (i.e., A1-A1 or A3-A3). Fowler et al. demonstrated this manner with an electron micrograph.¹⁵ Note that AA2MM discussed in this work focuses more on details within each of the clusters shown in the images from electron micrograph. AA2MM in this study is very similar to the prior model, except for the treatment of the A2 domain [see Fig. 2(a) as well as the text below]. Details of the model are presented below with emphasis on the different manner in which the unfolding and mechanical response of the A2 domain are described.

In this work, we performed Brownian dynamics simulations in an implicit solvent model with hydrodynamic interactions included.^32,33 In such simulations, the dynamics of the i-th bead position r_i is given by the following Langevin equation:³⁴ 

where $rit$ and $rit+Δt$ are the position of the i-th bead at time step t and t + Δt, respectively. In this work, each bead is characterized by bead radius a, and solvent flow is characterized as simple shear flow (or elongational flow) with a constant viscosity η (all parameter values used here can be found in Table I). Based on the assumption that solvent flow properties, i.e., primary velocity field will not change during the simulation process, all the solvent parameters, along with the bead radius, are represented by a drag coefficient ζ (see Table I). Additionally, it is customary to present the equations of motion in a dimensionless form where the length, force, and time are scaled by the bead radius a, k_BT/a, and ζa²/k_BT, respectively.^35,36 k_B is the Boltzmann constant, and T is the temperature. ∇u is the velocity gradient tensor that describes the primary flow field of the solvent. Polymers are immersed in fluid condition with the undisturbed velocity of the solvent felt by a given bead as v^∞ = (∇u)^T · r, where r_i is the position of bead i. For simple shear flow where the vorticity direction is along y, v_i can be calculated through $v∞=γ̇zx^$⁠, where $γ̇$ is the fluid shear rate in pure shear flow, _z is the coordinate in the z direction, and $x^$ is the unit vector in the x-direction.²⁸ Under pure elongational flow, $v∞=ϵ̇(zz^−12xx^−12yy^)$⁠, where $ϵ̇$ donates the strain rate, x and y are the coordinates in x and y directions, and $y^ and z^$ are the corresponding unit vectors.³⁷ Hydrodynamic interactions are manifested in D_ij, which represents the diffusion tensor and is calculated by

where δ_ij is the delta function, I is the identity matrix, and Ω_ij is the hydrodynamic interaction tensor. In this work, we employed the Rotne-Prager-Yamakawa (RPY) tensor for Ω_ij.^38–40 When simulations are performed with molecules in proximity to a wall (i.e., a surface positioned at z = 0), a modification to Ω_ij is employed to ensure the no-slip condition at the infinite planar wall.^41,42 This means hydrodynamic interactions explicitly depend on the bead heights in the z-direction for an infinite flat wall in the xy-plane. This modification term is identically zero in the absence of a wall.³⁶ $(∑j=1N∇rj⋅Dij)Δt$ is a displacement arising from the divergence of the diffusion tensor. It accounts for a spurious bead flux that manifest as a result of satisfying the no-slip condition at the planar wall. The divergence of the diffusion tensor is given by⁴¹ 

In Eq. (1), $6Δt∑j=1iλij⋅nj$ represents the Brownian random displacement, where n_j is a random three-dimensional vector and λ_ij is calculated by $δijI+ζ∑j=1NΩij=∑l=1λilλlj$⁠. Specifically, Cholesky decomposition was used to obtain the magnitude of λ_ij via $λii=Dii−∑p=1i−1λipλip$ when i = j and $λii=1λii(Dii−∑p=1i−1λipλjp)$ when i > j.³⁴ 

In Eq. (1), F_ij represents the force acting on i-th bead due to the j-th bead, which includes a bead-bead interaction force (⁠$FijLJ$⁠) that acts between all pairs of beads in the system. This force term helps limit bead-bead overlap; furthermore, the interaction energy parameters for this term are chosen so that multimers in quiescent and low shear rate flow conditions adopt a compact conformation. Adjacent beads in a multimer are also subject to one of two possible bead-bead bonded force terms, depending on whether the beads are in different monomers or the same monomer (⁠$Fijharm$ or $FijA1−−A3$⁠, respectively, see below),

The following Lennard-Jones (LJ) force is used for $FijLJ$⁠:

where ε and σ are the energy and length parameters. max(a, b) is the rectified linear unit function that returns the max value of a and b. In this work, the max function is used to plateau the repulsive Lennard-Jones force at very close separation, which corresponds to 80% overlap of the beads and smaller (i.e., 1.65 in reduced units). For nonbonded beads, $FijLJ$ is the only nonzero term in Eq. (3). For the additional force terms between bonded beads, A1-A1 and A3-A3 (i.e., intermonomer) connections are represented by relatively stiff harmonic springs with the same spring parameters,

where f_harm is the spring constant and $r̃harm$ is the equilibrium spring elongation. The connection between A1 and A3 beads in the same monomer via the A2 domain will be described in detail in the section titled Activated A2 domain model. Thus, in this work, A1/A3 beads are treated identically; however, when modeling tethering of a vWF multimer to a collagen-coated surface, only A3 beads are permitted to bind to designated collagen binding sites, in concert with what has been observed experimentally for vWF.

As described above, our prior model characterized the mechanical response of the A2 domain via a FENE spring. Furthermore, the force and length parameters of the FENE spring were based on the Worm-Like-Chain (WLC) force-extension correlations of unfolded A2, as obtained from experiment.^25,29 In our prior model, unfolding was not explicitly represented but was instead associated with FENE spring force excursions that were higher than a force threshold, as guided by experiments.²⁸ Even though this model is suitable for representing the extendable property of A2 domain, it does not explicitly represent the change (Σ^f→u) from a folded (Σ^f) to an unfolded (Σ^u) state or vice versa (Σ^u→f).^25,28,29 Chemically, the transition Σ^f→u or Σ^u→f is an activated process where the energy landscape is influenced by external forces (i.e., mechanoreactive). Furthermore, it is associated with a discontinuous force vs extension behavior that did not exist in our prior model. In order to explicitly include unfolding in our simulations, we introduce here a new activated A2 monomer model (AA2MM).

The goal in this work is to capture the mechanoreactive nature of the vWF A2 domain during Σ^f→u or Σ^u→f process; thus, we need to characterize in some way the influence of A2 domain distortion on the energy landscape. From force spectroscopy experiments, as shown by Zhang et al.,²⁵ unfolding was visualized by a sudden force drop with little change in the extension. The instantaneous force for Σ^f→u (i.e., F^f→u) obeyed a normal distribution for each force loading rate; furthermore, the most probable unfolding force was found to increase with rate. After unfolding, if a constant force was maintained on the domain, an additional extension of order 50 nm was observed. The authors also probed refolding kinetics for a range of applied force on the domain and found that refold probability for a given amount of time decreased with increasing tension. Zhang et al.²⁵ fit data from their pulling experiments to a rate law $kuF=ku0 exp(F/fβ)$⁠, where k_u is the unfolding rate, $ku0$ is the unfolding rate in the absence of applied force, F is the force on the domain, and f_β is a normalizing force parameter. Refolding data were well described by a similar model,²⁵ but with different $kf0$ and f_β. Note that, in the main manuscript of Ref. 22, refolding data were instead fit to a model by Evans that accounts for protein compliance; however, in the supplementary material for that work, the authors compared fitting refolding data with the Bell model to the Evans model and the differences were small. A similar comparison of these models in a related study was presented more recently.⁴³ Thus, here, the same model with different parameters is used for both unfolding and refolding. These rate laws are reflective of a Bell model⁴⁴ for binding in biological systems where applied external force can shift the energy landscape dictating the unbinding/binding process; here, the activated process is unfolding/folding, but the effect of force on the energy landscape is well described by a Bell model. For Σ^f→u or Σ^u→f of A2, Fig. 1 gives a schematic diagram of relative free energy levels between folded (Σ^f) and unfolded (Σ^u) states depending on external forces. If no external force is applied, Σ^f has a lower energy level compared to Σ^u, and the corresponding energy landscape is shown with a blue solid line in Fig. 1. Applied external force changes the relative free energy level of Σ^f and Σ^u, making the A2 domain easier to unfold. If the applied external force is large enough, the unfolded state can become energetically favorable to the folded one; such a scenario is depicted by the orange dashed line in Fig. 1.

To incorporate mechanically mediated unfolding/refolding in our CG model, we first define A2 domain distortion based on the edge-to-edge separation of the A1 and A3 beads in a monomer $DEE=rij−r̃A1−−A3$⁠, where $r̃A1−−A3=1$ nm represents the equilibrium edge separation for adjacent A1 and A3 beads. Prior work by our group as well as others explored A2 domain behavior using all-atom molecular dynamics simulations, and the undistorted morphology of the A2 domain was found to be very approximately spherical with a diameter ranging from 3 nm to 4 nm^26,27. However, it is known from experiments that the vWF monomer structure exhibits a smaller separation distance between the A1 and A3 domains because the majority of the A2 domain does not locate directly between those domains. The equilibrium $r̃A1−−A3=1$ nm used here reflects the equilibrium separation distance between A1 and A3. A positive D_EE represents a tensile distortion of the A2 domain, while negative means compressive distortion although negative excursions are rare and of small magnitude due to the $FijLJ$ acting between beads. In our AA2MM, a system of two springs is used to describe the mechanical response of the A2 domain [Fig. 2(a)]: a FENE spring (S_A for convenience) with the same parameters used in our prior monomer model and a relatively stiff harmonic spring (S_B). In the folded state, both springs act in parallel; however, the harmonic spring does not engage unless D_EE exceeds a designated slip length of A2 domain (L_slip) that is set to be 2 nm in this work [Fig. 2(b), blue shaded region]. This is to model the relatively loose force-extension correlation at a very small separation distance of the cluster of A1, A2, and A3 domains within each monomer; such behavior can be associated with domain rearrangement that occurs prior to force being applied on specific chemical bonds. Therefore, when D_EE is smaller than L_slip, the harmonic spring will not be engaged and the force-extension correlation at very small tensile distortion is dictated by the FENE spring [Fig. 2(b), blue shaded region]. When D_EE becomes larger than L_slip, the harmonic spring engages [Fig. 2(b), orange shaded region] and, due to its relatively stiff nature, dominates force vs extension until unfolding occurs. Therefore, $FijA1−−A3$ is

where f_A and r_max are the FENE spring constants and f_B is the harmonic spring constant. The max function permits implementation of the slip length criterion. $Πu$ is a Boolean operator; for each A2 domain in the system, $Πu$ = 0/1 indicates the domain is folded/unfolded, respectively. The algorithm for determining whether or not a state change occurs is described below. If an A2 domain is unfolded (⁠$Πu$ = 1), the harmonic spring is considered broken such that only the FENE spring acts between adjacent A1 and A3 beads.

Figure 2(c) shows the Σ^f→u decision tree in AA2MM. For each A2 domain in a system, when D_EE is smaller than L_slip, the A2 domain is considered to be in a folded state and, as described above, the harmonic spring is not engaged. If a given A2 domain is in the folded state and D_EE > L_slip, a force dependent probabilistic check is made to determine if that A2 domain undergoes Σ^f→u or Σ^u→f. If the check is successful, on the following time step, that A2 domain is considered unfolded [i.e., $Πu$ in Eq. (6) above is changed from 0 to 1]. To determine the force dependent probability, the same rate law that was used to fit experimental data is used here. For unfolding, $kuF=ku0 exp(F/fβ)$⁠, where F is given by the instantaneous value of $FijA1−−A3$⁠; the same expression, with different parameterization, is used for folding. The parameters for unfolding, $ku0$ and $fβu$⁠, and for folding, $kf0$ and $fβ f$⁠, are taken directly from experiment.²⁵ To convert the rate to probability, k_u(F) [or k_f(F)] is divided by the attempt frequency. Thus, when an attempted state change is modeled, a random number uniformly distributed between 0 and 1 is selected and, if the number is less than the computed probability, the state change occurs and $Πu$ is switched from 0 to 1 (unfolding) or from 1 to 0 (folding) in Eq. (7). The simulation time step (⁠$δt$ = 0.1 μs) ensures that D_EE varies continuously; nonetheless, the sudden disengagement of the harmonic spring that occurs upon unfolding produces a discontinuous decrease in $FijA1−−A3$⁠, which is schematized with a blue dashed line in Fig. 2(b). This is analogous to the force drop upon unfolding observed experimentally. For unfolded A2 domains in the model, if D_EE between its two adjacent A1 and A3 beads drops below L_slip, the A2 domain is considered folded for the subsequent time step.

All parameter values used in this study are summarized in Table I; in all cases, parameter values were either taken directly from experimental data, known structural characteristics of vWF or taken from insight obtained from fully atomistic simulations. The radius for both A1 and A3 beads (a in the expression for drag coefficient) is based on the known structure of the vWF monomer. ε and σ in Eq. (4) were set to give an equilibrium D_EE of 1 nm between adjacent A1 and A3 beads in a monomer based on experimental structural data. Furthermore, the chosen ε value causes multimers in low shear rate flow conditions to adopt a compact conformation. The length parameter in the harmonic bonds connecting monomers [$r̃harm$ in Eq. (5)] was also chosen based on the known structure of vWF, and the force parameter [f_harm in Eq. (5)] is chosen to be somewhat arbitrarily large to model the strong intermonomer bonds in vWF. All parameters in Eq. (6), except for L_slip, were taken directly from Ref. 29. The value L_slip = 2 nm was selected somewhat empirically by considering the A2 domain length distribution due to random bead displacements in quiescent flow simulations with no slip length and all A2 domains always in Σ^f (Fig. S1A). The value is small compared to the bead radius in our model; furthermore, our prior atomic scale simulations of A2 domain pulling showed no unfolded structures within this L_slip (Fig. S1B). Sensitivity analysis of f_B and L_slip will be discussed later in detail. All parameters in both the unfolding and folding rate equations were again taken directly from Zhang et al.²⁵ It has previously been demonstrated that the attempt frequency for state changes in complex macromolecular systems such as proteins is orders of magnitude lower than the attempt frequency for activated processes involving a small number of atoms. For protein unfolding, attempt frequency of order 10⁷ Hz has been suggested.⁴⁵ Since the time step in this work is 10⁻⁷ s, a state change attempt was made every 100 time steps; however, because the rate of state change is normalized by the attempt frequency to obtain the probability, results presented herein do not change if attempts are instead made every 10 or 1 time steps.

To validate our AA2MM implementation, we conducted single monomer pulling simulations wherein, for each condition tested (e.g., each loading rate), 2000 independent monomers were simulated. The external force acting on the A2 domain in each case was increased at a constant rate set by the loading rate. For each A2 domain modeled, the value of F^A1–A3 [Eq. (7)] at which the state change Σ^f→u occurred was recorded as the unfolding force for that domain. Separate monomer simulations were performed (2000 independent simulations in each case) to assess refolding statistics. Specifically, unfolded A2 domains were held at a constant force, and the fraction of refolded domains was recorded as a function of time.

Results from monomer simulations are shown in Fig. 3. The most probable unfolding force vs loading rate is shown in Fig. 3(a), and a very satisfactory agreement with experiment can be seen. Figure 3(b) shows the probability distributions of the unfolding force for each loading rate modeled. These can be compared to the histograms for unfolding force presented in Ref. 25, and a very satisfactory agreement is obtained. To illustrate such a comparison, the inset of Fig. 3(b) shows the probability distribution for an intermediate pulling rate (35 pN/s) plotted similarly as data presented in Ref. 25, along with corresponding data in Ref. 22. It can be seen that the corresponding experimental data show a little more spread in unfolding force; however, in simulations, rather than deriving force on the domain indirectly from the optical tweezers, the force is directly computed from the coarse-grain molecular ensemble. Thus, the systematic variation in the measured force exists in experiments but does not in the simulation. Finally, Fig. 3(c) shows results from simulations for unfolded domains subject to constant holding force. For the two lower values of holding force, an agreement between experiments and simulations is again highly encouraging. At the highest holding force, there is only a single data point from experiment; although an agreement between simulations and experiment is not as satisfactory for those data, it can be seen that the general trends in these experimental data are also well captured by AA2MM.

Having compared predictions from the new monomer model with experimental data that inspired the model (and were used to parameterize the model), we next examined vWF multimer behavior as predicted by AA2MM. vWF multimers were modeled in three flow scenarios: tethered to model collagen-coated surfaces while subject to shear flow, in bulk shear flow (i.e., far from any surface), and in bulk elongational flow. For all scenarios, the influence of the shear/elongation rate on multimer behavior, including submonomer unfolding behavior, was studied. For all multimer simulations, 10 independent single-multimer simulations were performed, giving 10 distinct samples for statistical analysis. For each simulation, the starting state of the polymer was formed by equilibration dynamics simulations from the same starting conformation with the first bead tethered to the collagen wall but conducted with a different random seed. The initialization for bulk flow simulations (described in the section titled Bulk flow simulations) used a similar procedure as the tethered ones (but with no collagen-coated wall). After this initial equilibration, flow simulations began at the desired shear or elongation rate.

The results of tethered molecule simulations are presented first to enable a comparison with recently published experimental data by Fu et al.³⁰ Here, we modeled 112-monomer chains with the first A3 bead in each chain tethered to the model of the collagen-coated surface. Polymer extension values under various shear rates were calculated for comparison with the experimental data from Fu et al.³⁰ All simulations began with polymers in a globular conformation (Fig. S1 C) before the flow was turned on. After shear flow began, polymers responded by gradually unraveling and extending along the flow direction until a shear rate dependent plateau in extension was reached; the plateau value increased with increasing shear rate (Fig. S2). In addition, with increasing shear rate, some A2 domains were observed to unfold. In this work, we computed the normalized chain extension in the flow direction to represent the polymer conformation. The normalized extension is defined as the measured length minus the initial length, divided by the difference of the maximum and initial lengths. Initial and maximum lengths are given by simulations at shear rates zero and 150 000 s⁻¹, respectively.

Figure 4(a) shows the average normalized extension values from simulations along with the experimental results by Fu et al.,³⁰ who utilized the same calculation. Results from the AA2MM model show very good agreement with the experimental data.³⁰ It should be noted that multimer elongation is a complex process influenced by many factors, i.e., A2 domain unfolding, intramolecular interaction between different domains, and also the conformational change in any part of the multimer. The extension discussed in Fig. 4(a) only focus on the multimer level irrespective of extension contributed from any specific factor. By applying AA2MM, a good match can be observed with experimental data at the multimer level. In addition, AA2MM has the ability to represent A2 (un)folding. Because all domains other than A2 are coarse-grained into one of the two beads making up each monomer in AA2MM, the influence of all domains other than A2 on multimer elongation is therefore represented in the bead-bead interaction. As such, this model cannot quantify the fraction of multimer extension due to distortion of submonomer entities other than the A2 domains. The specific roles of other domains in controlling the mechanical response of vWF need further investigations. A2 domain pulling experiments using optical tweezers⁴⁶ provide some insight into the appropriate stiffness of response prior to A2 domain unfolding; however, variation is observed. From those data, reasonable bounds on the stiffness parameter are 1 pN/nm ≤ f_B ≤ 5 pN/nm and this stiffness parameter has a little influence on results shown here in Fig. 4(a). Nonetheless, given the CG nature of the model, agreement with experimental data is very satisfactory.

To explore the A2 unfolding behavior in tethered multimers subject to shear flow, we defined a polymer unfolding probability. To be clear, this does not refer to the overall multimer conformation (i.e., globular vs partially unraveled vs fully or mostly unraveled). Here, “polymer unfolding” is used differently: that is, any single multimer is considered unfolded if at least one of its A2 domains is unfolded; given this, the percentage of time for which this was true gives the polymer unfolding probability. Again, simulation results were averaged over 10 statistically distinct samples. From Fig. 4(b), it is clear that negligible unfolding occurs at zero flow. However, unfolding probability increases to approximately 60% for $γ̇=2500s−1$ and it increases sharply again to ∼90% at $γ̇=5000s−1$⁠. Polymer unfolding probability reaches almost 100% when the shear rate equals 25 000 s⁻¹.

From tethered multimer simulations, the dependence of domain unfolding probability on the A2 position along the multimer was explored. Similar to polymer unfolding probability, domain unfolding probability was computed as the percentage of time each domain along a multimer was unfolded. Figure 4(c) shows domain unfolding probability as a function of the A2 position along the polymer contour. It is clear that A2 domains located near the tether point are more likely to unfold, which agrees with the notion that hydrodynamic forces are greatest in that region. When $γ̇=2500s−1$⁠, domain unfolding probability is ∼60% near the tether point and then decreases to zero for domains that are more than ∼60 monomers from the tether point. For $γ̇=5000s−1$ and higher shear rates, two distinct regimes of decreasing probability are observed. For example, for $γ̇=5000s−1$⁠, the domain unfolding probability near the tether point is ∼90%. This again decreases with increasing distance from the tether; however, it initially decreases more gradually than it does for monomers that are more than ∼70 from the tether point. Put differently, a “knee” exists in the $γ̇=5000s−1$ curve at ∼70 monomers from the tether point; no unfolding occurs for monomers more than 90 away from the tether point. For higher shear rates, the position of the knee in the curve increases in both monomer position and unfolding probability. For $γ̇=25000s−1$ and higher, a very high probability for domain unfolding persists essentially all along the chain. In the absence of unfolding, the distribution of force experienced by a tethered molecule exhibits a fairly consistent rate of decrease as the distance from the tether point increases; this is in contrast with behavior observed here. We associate this with the stress relief mechanism imparted by domain unfolding. That is, unfolded domains are much more compliant than folded ones; thus, unfolding along the part of the chain closer to the tether point reduces force experienced by segments of the multimer farther from the tether. Domain unfolding behavior shown in Fig. 4(c) bears out that this results in a more abrupt change in the way force varies along the tethered multimer, particularly beyond the knee in the curve.

Although data in Fig. 4(c) cannot be compared directly to data from Fu et al., a qualitative comparison can be made between our domain unfolding probability trends and observations from their experiments on binding to tethered vWF multimers. Specifically, in that paper, GPIbα was introduced into the system, while multimers were subject to shear flow. This entity is known to bind to the A1 domain of vWF; because they are colored with different fluorescence signals from vWF multimers, their fluorescent signal provides spatial information on where binding is occurring. Most notably, multiple groups have advanced data supporting the notion that such A1- GPIbα binding is facilitated by the A2 domain adjacent to A1.³⁰ That is, binding is only possible if the A2 domain is subject to some degree of tensile force, perhaps sufficient to drive unfolding but the amount of distortion required is unclear. From experiments by Fu et al.,³⁰ images of bound fluorescent GPIbα molecules showed a preference for binding along the part of the chain closest to the tether point with the end of the chains farthest from the tether typically showing no or little binding. This observation at least qualitatively agrees with the trends in domain unfolding probability obtained from simulations here. A separate qualitative comparison can be made regarding A2 domain-ADAMTS13 blood enzyme binding, which leads to proteolysis of vWF molecules and serves as a critical size maintenance mechanism for vWF in blood. It has been stated in a collection of previous experimental works^25,30,46 that such binding should only occur for shear rate $γ̇=5000s−1$ and higher, and A2 unfolding precedes such binding; trends observed here in domain unfolding with varying shear rate agree at least qualitatively with that assertion.

Simulations in bulk flow conditions were performed for pure shear flow and pure elongational flow with polymers consisting of 80-monomers. We calculated a radius of gyration (Rg) to represent the conformation of each polymer under varying shear/strain rate bulk flow conditions,

where r_i represents the position of each bead, while r_com is the center of mass position of the whole polymer. Angle brackets denote an ensemble average over N = 160 beads in each chain, for all chains modeled, for all time. Average Rg vs shear rate is shown in Fig 5(a).

From Fig. 5(a), it is clear that Rg increases with the shear rate similar to what has been observed for previous CG models of vWF, including our group’s prior model.^28,29,36 Indeed, significant simulation data exist characterizing the behavior of CG vWF multimers in shear flow and results with the new model presented here generally bear out the same conclusions. For example, for globular polymers, a protrusion-nucleation theory has been advanced where the mechanism for polymer unraveling is the protrusion of a section of the multimer away from the compact globule. The length of the protrusion and its orientation to the flow field then dictate whether further unraveling occurs or if the protrusion instead retracts into the globule.²⁸ Similar behavior is exhibited by polymers modeled here; to demonstrate this, conformations are shown in supplementary material, Fig. S1 C. It has also been observed in prior simulations that, for high shear rate where complete or nearly complete polymer unraveling occurs, multimers sometimes exhibit rapid dynamic transitions between more globule conformations and unraveled conformations. Again, the model presented here exhibits similar behavior.

Figure 5(a) reveals significantly different behavior in elongational flow than what is observed in shear flow. The transition from compact to elongated is much more abrupt in elongational flow, and it occurs at a much lower elongation rate. Indeed, the approximate elongation rate demarking the transition from compact to elongated conformations is $ϵ̇=3500s−1$⁠, which is reasonable with regard to physiological blood flow conditions. That is, the multimer size studied here would not be expected to elongate in healthy blood flow conditions and the transition shear rate identified is high with respect to normal flow but might be realized near vascular injury locations. Elongational flow is significantly more effective in causing molecular elongation than pure shear flow because multimer tumbling, which is a common behavior exhibited by compact multimers in pure shear flow, is suppressed in elongational flow. This has been observed by prior authors.³⁷ Results here provide further evidence that, for a multimer within the observed healthy size distribution for vWF, elongation would not be expected in either shear or elongational flow in typical blood flow conditions. In regions of relatively high elongation rate (i.e., near injury), elongation is predicted.

The monomer model advanced here is distinct from prior work in that we can explicitly identify when in the model unfolding occurs. Similar to the analysis done for tethered molecules, we computed the polymer unfolding probability for multimers in bulk flow. Data in Fig. 5(b) show that this probability remains close to zero until the shear rate exceeds $γ̇=200000s−1$⁠. For $γ̇>350000s−1$⁠, the unfolding probability increases more significantly such that at the highest shear rate explored, the unfolding probability is ∼70%. It is interesting to note that when the shear rate reaches 350 000 s⁻¹ or higher, the predominant dynamic behavior of multimers changes. That is, instead of the tumbling of globular/protruded conformations, polymers start to show some fully unraveled conformations (Fig. S1 C), which exposed more A2 domains to the hydrodynamic force from shear. Associated with this conformational behavior transition is a transition in the distribution of force on A2 domains as a function of position along the multimer. Specifically, globular conformations are associated with a peak in force near each end of the polymer and the peak force value is relatively low; unraveled multimers exhibit a single peak in the force profile near the center of the multimer chain and at a peak tension that is five to six times the peak force for globular conformations. These observations with the new model presented here are very similar to what was presented in a prior work by some of the current authors.²⁸ Nonetheless, force distribution data bearing out this discussion for the model argued here are shown in Fig. S3 of the supplementary material. This dynamic conformational (and force distribution) transition is what leads to the more sharply increasing unfolding probability. In elongational flow, for elongation rates higher than $ϵ̇=3500s−1$⁠, unfolding transitions from being a very rare event below this elongation rate threshold to essentially 100% unfolding probability once multimer unraveling is predicted.

Parameter sensitivity, with regard to the A2 domain spring parameters that are newly introduced here, has also been tested for pure shear flow simulations. For f_B, values of either 1 pN/nm or 5 pN/nm were used; for L_slip, values of either 2 nm or 3 nm were used. For all four combinations of these two parameters, a very little difference was observed for Rg as a function of shear rate, indicating that Rg is more dependent on the Lennard-Jones parameter ε than A2 domain spring parameters chosen. For polymer unfold probability (Fig. 6), which represents monomer scale behavior, there is a little change for f_B = 5 pN/nm regardless of the value of L_slip. However, reducing f_B lowers unfold probability; this is directly related to the dynamics of D_EE. When an A2 domain is in the folded state, periodic excursions of D_EE larger than L_slip sample the harmonic component of the A2 domain spring system; lowering f_B reduces the force fluctuations associated with D_EE fluctuations, thus lowering the force-based probability for unfolding to occur. Perhaps nonintuitively, the probability is reduced for lower L_slip. We attribute this to a momentum effect. When two beads experience a close approach, the bead-bead interaction generates repulsion between them, causing them to move apart; deceleration due to bead-bead attraction does not become significant until D_EE exceeds L_slip. Thus, larger L_slip permits the beads to achieve a greater velocity before significant deceleration manifests, permitting them to reach larger separation distance, with concomitantly larger tensile domain force. Such excursions lead to a greater unfold probability at a given shear rate for larger L_slip.

To further investigate the unfolding event of the A2 domain, we calculated the duration of single A2 domain stays in Σ^u; each unfolding event was counted during a total period of 15 × 10⁻⁵ s. Then, these unfolding events were separated into 10 bins with the time interval of 15 × 10⁻⁶ s [Fig. 7(a)]. From this figure, it is clear that the first bin (≤15 × 10⁻⁶ s) dominates the count and only the first 4 bins contain nonzero values. For each bin [Fig. 7(b)], the count profile was shown as a function of the A2 position along the contour, which yields a different profile shape for the first bin compared to the rest. For the first bin, which represents the shortest unfolding durations, a double peak curve for count vs A2 domain position results with peaks near the multimer ends. On the other hand, when unfolding duration is higher than 15 × 10⁻⁶ s, the distributions show higher counts for A2 domains more central in the polymer chains. These results correlate with the force distributions vs A2 domain position discussed above (and shown in Fig. S3). That is, the longest duration unfolding events are associated with those domains on which the largest tensile forces are observed to act (i.e., for unraveled multimers and nearer to their center). Regarding the ADAMTS13 size regulation process of vWF polymers, even though the kinetic process has not been fully revealed yet, data in Fig. 7(b) may elucidate some further insight into the ADAMTS13-vWF activation time scale. For example, it has been showed that the ADAMTS13-vWF scission process may be associated with a larger time scale compared with the A2 unfolding or A1-platelet binding processes.^47,48 Thus, it is reasonable to believe that these different profile shapes in Fig. 7(b) may have some relation to the ADAMTS13-vWF size regulation process. Specifically, the relatively longer timed ADAMTS13-vWF scission preferentially occurs in the central region of unraveled vWF proteins.

Simulation predictions using AA2MM exhibit good agreement with experimental data from both monomers pulling experiments and microfluidic chamber flow-based experiments on vWF multimers. Simulations in various bulk flow conditions predict that vWF is most likely driven into an elongated state via elongational flow. Furthermore, for the multimer size studied here, which is within the expected size range for vWF in healthy human conditions, elongation is only predicted for shear flows that are elevated relative to healthy flow conditions. Such conditions may occur in the vicinity of vascular injury, which is where vWF elongation would be needed to abet the clotting cascade by binding to both collagen on vessel walls and platelets in the bloodstream. The new model, while very coarse-grained, thus seems to accurately characterize the mechanoreactive nature of vWF on both the multimer and submonomer scale.

It is of interest to further consider mechanisms of size regulation of vWF by the blood enzyme ADAMTS13. It is known that vWF molecules, when initially secreted into the blood, are ultralarge molecules that are many times the size of functional vWF molecules in blood.¹³ Depending on the method of assessment, a very approximate range of functional vWF size in healthy blood is from 5 to 200 monomers; some sources place the upper limit at 100 monomers or less.³ Regardless, the molecular weight of vWF proteins when they first enter the blood is significantly larger than indicated by either of the quoted upper bounds. If vWF remains at such large sizes, it will cause thrombosis by binding platelets in normal blood flow conditions; this is what occurs in some vWF-related diseases such as thrombotic thrombocytopenic purpura. Thus, quickly after secretion, their size is reduced by scission, caused by ADAMTS13 binding to unfolded A2 domains. That is, the binding site for ADAMTS13 within the A2 domain is only exposed after unfolding, and domain unfolding has therefore been advanced as the rate-limiting step for the scission reaction. It is further believed that multimer elongation precedes A2 domain unfolding because this permits hydrodynamic forces to more directly act on the A2 domains in an elongated vWF protein. Ultralarge vWF initially secreted into blood experience significantly greater forces due to hydrodynamic effects: they are more prone to multimer unraveling as well as submonomer domain unfolding. Thus, soon after entering the bloodstream, ultralarge vWF become elongated, A2 domains unfold, and scission occurs. From fragments so formed, if they remain larger than the functional size range, they too will relatively quickly experience multimer elongation, A2 domain unfolding, and scission. In this way, the combination of healthy blood flow conditions and the mechanics of vWF response—both at the multimer scale and at the submonomer scale—ensure that an appropriately sized population of vWF is established within a relatively short time (minutes to hours) after secretion into the bloodstream.³ Even for vWF multimers in the healthy size range, the lifetime of a vWF protein is limited, and it is believed that the recycling of normal-size vWF multimers also occurs via ADAMTS13 catalyzed scission albeit at much slower rate. Prior authors have combined coarse grained simulation with experiment to explore the flow activated scission by ADAMTS13⁴⁹.

Simulation results presented herein for pure bulk shear flow show negligible A2 domain unfolding until very high shear rates, much higher than physiologically relevant ones. This is partly due to a very low percentage of vWF molecules undergoing significant elongation until very high shear rates. Similar to what has been shown by prior authors,³⁷ dramatically different behaviors are exhibited by model vWF in elongational flow conditions. Elongation flow is experienced in the bloodstream near narrowing of the vasculature; however, because blood is a complex fluid, the flow should most likely be considered a combination of shear and elongation where the amount of each varies spatially and temporally. This type of flow is much more effective in driving multimer elongation, and this is convincingly shown by data presented here. Tumbling is a primary dynamic mode of vWF multimers observed in shear flow, but elongational flow effectively prevents tumbling, causing vWF multimers to exhibit a higher energy response mechanism (i.e., elongation) at lower rates of shear. Results from simulations here provide convincing evidence that vWF elongation (and subsequent A2 domain unfolding) preferentially occurs in elongational flow regimes. As observed above, domain unfolding is not observed here in low shear rate conditions akin to healthy blood flow—for either pure shear or elongational flow; however, this is expected for the size multimer modeled. Most notably, if scission requires relatively long-lived unfolded A2 domains, the data presented here provide strong evidence that such scission occurs nearer to the central region of elongated proteins. Such a self-similar scission process has previously been suggested as a mechanism to explain the banded nature of vWF size assays obtained, for example, from gel filtration. Results presented here dovetail well with this argument.

In this work, we present a new coarse-grained model for the vWF monomer, AA2MM, to represent the probabilistic unfolded/refolded state change in A2 domains in vWF polymers, as well as the change in mechanical properties of a monomer associated with changes between the unfolded and folded states. A system of springs is employed to represent the A2 domain force response under extension, and a probabilistic process, whose parameters are taken directly from domain pulling experiments, is used to determine when domain unfolding (or refolding) occurs. AA2MM has been validated through the comparison of monomer and multimer scale simulations with experimental results.^25,30 Results from AA2MM show good agreement with data from experiments. Finally, we implement AA2MM to predict the vWF polymer behavior and unfolding properties in bulk flow. Bulk flow A2 unfolding is intrinsic to the size regulation process of ultralarge vWF that are initially secreted into the blood; it is also related to some vWF-related diseases. Data presented here give further evidence that long-lived unfolded A2 domains preferentially occur nearer to the center of unraveled vWF proteins. This supports the argument that the banded nature of vWF size distributions is a result of scission via blood enzyme ADAMTS13 interaction with unfolded A2 domains preferentially occurring nearer to the center of the protein chains.

In the supplementary material, Fig. S1 shows A2 length distribution in quiescent flow. Moreover, Fig. S1 presents the snapshot of A2 domain structure and different polymer conformations. Figure S2 displays chain extension length as a function of time for tethered simulation. Averaged A2 domain force profiles are described in Fig. S3.

This work was supported by the National Science Foundation, Grant Nos. DMS-1463234 and DBI-1707207. Compute resources were provided by the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by the National Science Foundation, Grant No. ACI-1548562.⁵⁰ Finally, support from the Lehigh sol cluster is also acknowledged.