Molecular Mechanism of Protein Unfolding under Shear: A Lattice

Jan 12, 2018 - Proteins are marginally stable soft-matter entities that can be disrupted using a variety of perturbative stresses, including thermal, ...
0 downloads 8 Views 4MB Size
Subscriber access provided by Gothenburg University Library

Article

Molecular Mechanism of Protein Unfolding Under Shear: A Lattice Boltzmann Molecular Dynamics Study Fabio Sterpone, Philippe Derreumaux, and Simone Melchionna J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.7b10796 • Publication Date (Web): 12 Jan 2018 Downloaded from http://pubs.acs.org on January 12, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry B is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Molecular Mechanism of Protein Unfolding Under Shear: a Lattice Boltzmann Molecular Dynamics Study Fabio Sterpone,∗,† Philippe Derreumaux,† and Simone Melchionna‡ †Laboratoire de Biochimie Th´eorique, IBPC, CNRS UPR9080, Univ. Paris Diderot, Sorbonne Paris Cit´e, 13 rue Pierre et Marie Curie, 75005, Paris, France ‡ISC-CNR, via dei Taurini, 00185, Rome, Italy E-mail: [email protected]

1 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 23

Abstract Proteins are marginally stable soft-matter entities that can be disrupted using a variety of perturbative stresses, including thermal, chemicals or mechanical ones. Fluid under extreme flow conditions is a possible route to probe the weakness of biomolecules and collect information on the molecular cohesive interactions that secure their stability. Moreover, in many cases, physiological flow triggers the functional response of specialised proteins as occurring in blood coagulation or cell adhesion. In this paper, we deploy the Lattice Boltzmann Molecular Dynamics technique based on the coarse grained model for protein OPEP, to study the mechanism of protein unfolding under Couette flow. Our simulations provide a clear view on how structural elements of the proteins are affected by shear, and for the simple study-case, the β-hairpin, we exploited the analogy with pulling experiments to quantify the mechanical forces acting on the protein under shear.

Introduction Important biological processes occur under fluid flow conditions and depend on the response of proteins to the mechanical stress arising from the surrounding fluid, causing shearing, elongational and tensile forces. A very pertinent example is the coagulation process triggered by a sudden increase in blood flow caused by a vessel injury. In homeostasis, the functional response of the von Willibrand factor (vWf) is associated to its elongation and the consequent exposure of the binding site to platelets activation 1 , and to the mechanical unfolding of one its constituting domain, the A2 domain, that allows the cleavage of the chain 2 . A second example is the action of catch-bonds proteins 3 that favour cell adhesions to target membranes in the vascular system or in the urinary tract. These proteins are activated by shear flow that induces allosteric-like conformational changes 4 . In recent years it has been also debated whether shearing forces can cause protein unfolding 5–9 . In fact, beside the activation of protein functionality described above, shear flow can 2 ACS Paragon Plus Environment

Page 3 of 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

perturb proteins in many technological processing (e.g. by filtering or pumping devices) compromising their activity. Moreover, it is now clear that shear flow promotes aggregation of misfolded and disordered protein systems 10–12 , although contradictory results were reported in the literature. Jaspe and Hagen 5 studied experimentally a small globular protein, the horse cytochrome C (104 amino-acids, a.a.), and proposed a simple but effective analytical model such that unfolding of small globular proteins is expected to occur only at very high shear rate, γ˙ ∼ 107 s−1 . In line with this prediction, it was recently reported 13 that the rhGH (∼ 200 a.a.) and IgGI (∼ 400 a.a.) proteins do not unfold up to highest tested shear rate of 104 s−1 . However, other investigations reported that important conformational changes can be observed even at low shear rate for a mid size protein, like lysozyme 6 (∼ 130 a.a. , unfolding observed at γ˙ ∼ 101 − 102 s−1 ), or a larger protein like Bovine Serum Albumin (BSA) 6 (∼ 580 a.a., unfolding observed at γ˙ ∼ 102 s−1 ). Helix to coil transition was also observed for the model peptide poly-L-lysine that loses its α-helical structure at very low shear rate 7,8 . Probably, the difficulty to have a complete microscopic and temporal control of the shear perturbation in experiments lays at the origin of such conflicting values for the critical shear rate. Computer simulation is a powerful tool to investigate the behaviour of proteins in a broad range of unfolding conditions (chemical, thermal and mechanical unfolding), and can provide atomistic insight in fluid flows as well. In order to explore the shear perturbation, various protein models have been used to obtain insights on the mechanism of unfolding under shear. In the specific case of the vWf system, a comparison with experiment, supported by scaling arguments, has been successful, and allowed to individuate the critical shear rate that triggers the transition from a compact globular configuration to an extended one 14,15 , with a vWf represented by a single polymer with attractive spherical beads. In other studies, 16–18 the protein is reduced to a sequence of carbon α atoms interacting with a G¯ o-like potential to favour the native configuration over non-native configurations. While the G¯ o model 19 can describe successfully some protein folding features, this model by not allowing favourable non-native contacts and

3 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 23

the side-chain specificity can miss important folding states 20,21 and therefore unfolding states under any conditions, and in particular under shear flows. The use of all-atom models in explicit solvent is instead rare, and very few studies attempted to monitor the effect of shear flow on local conformational changes in proteins involved in the blood coagulation process 22,23 . In this work we present an extended, and explorative, study based on the high-resolution coarse-grained model OPEP (Optimized Potential for Efficient peptide structure Prediction) 24 and the Lattice Boltzmann Molecular Dynamics technique 25 . The novelty of this approach, as compared to previous ones, is that the OPEP force field has atomistic resolution for the backbone and allows to describe protein folding and unfolding in high detail. Moreover, the Lattice Boltzmann technique naturally allows to generate the shear flow instantaneously coupled to the protein conformational dynamics, as we have tested in the preliminary investigations of the FimH and A2 vWf proteins under Couette flow 26,27 . Since the OPEP force field is also optimised to reproduce thermal unfolding in semi-quantitative agreement with experiments 28–30 , our simulations of shear-induced unfolding provide important indications on both the mechanism of unfolding and the tensile forces perturbing the protein native state. In fact, for a set of model proteins of small size we were able to observe and follow unfolding in the nanosecond timescale at very high shear rates (> 109 s−1 ) by considering multiple independent events. Moreover, for the simplest system studied here (the model β-hairpin) we explicitly link the unzipping mechanism to the drag forces generated by the shear gradient and periodically acting on the folding reaction coordinate as as result of the rotational component of the shear field. The obtained results encourage the application of LBMD simulations based on the OPEP CG model to other systems for which experimentally the shear induced unfolding has been reported, or that are functional under shear flow, e.g. the vWf domains.

4 ACS Paragon Plus Environment

Page 5 of 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Methods Simulations Lattice Boltzmann Molecular Dynamics (LBMD) simulations, used in the past to study colloidal and polymeric systems 31–33 , are here used, in combination with the OPEP force field 25 , to study protein unfolding under shear flow. In LBMD, a molecular model based on a implicit solvent force field is coupled to a kinetic representation of the solvent, simulated via the Lattice Boltzmann (LB) technique 34 . This technique enables the inclusion of hydrodynamic interactions in the context of coarse-grained Molecular Dynamics and to follow the local structure of the flow field in time. The coupling between molecules and the solvent has the form of a Stokes-like drag force acting on a per-particle basis:

F~iD = −Γ(~vi − u˜i )

(1)

where ~vi is the velocity of the particle i, ~u˜i is the fluid velocity ~u smeared over a finite extension of the i-th particle (see 25 for details). Γ is the friction, an adjustable parameter in the methodology taken to be equal for all particles. The drag force adds up to the usual conservative forces derived from the Hamiltonian of the system, FiC = −∇i U ({r}), and a random white noise, FiR , represents thermal fluctuations. The reader can find the technical details in Ref. 25 . For all simulations, we chose Γ = 0.1 f s−1 to ensure temporal decoupling between the timescales relative to the shear and the frictional forces. The LB evolution was solved by using the BGK (Bhatnagar-Gross-Krook) collisional operator 35 , with a lattice grid spacing of 2 ˚ A, sufficient to resolve local hydrodynamic interactions for macromolecular systems, the interested reader can find in our previous works 25,36 a detailed discussion about the size of the grid spacing used to couple LB to a quasi-atomistic bio-molecular force fields. The molecular and LB dynamics were evolved synchronously using a physical timestep of 1.5 f s. In order to

5 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 23

Z β-hairpin

WW domain

X

Ca domain

Figure 1: Left side: schematic representation of the Couette flow generating a velocity gradient along the Z direction, and an embedded protein (green body) in the simulation box. Right side: secondary structure cartoon of the three proteins considered in this work. generate the Couette flow with a velocity gradient oriented along the Z direction, a few layers of fluid at the opposite edges of the simulations box in the Z direction were uniformly accelerated using a constant directional force (f X ) along -X and X, respectively (see Figure 1). As a result of this drag, a fluid velocity gradient is naturally established along the Z axis. Using different magnitudes for the force f X , several shear rates γ˙ were scanned in the interval 109 − 1011 s−1 . The proteins were placed in the middle of a simulation cubic box of size L = 50 ˚ A and equilibrated in absence of shear. From the equilibration phase independent configurations with different orientations were selected for the simulations. In order to confine the system in the Z direction, and therefore avoid an artificial velocity discontinuity for a protein moving across the simulation cell boundaries, the edges of the box at Lz = 0 and Lz = L were represented via a static repulsive wall where the no-slip boundary conditions were imposed. It is worth remarking that, while in the generally used simplified scheme a static velocity gradient is imposed onto the protein, in LBMD the velocity gradient responds dynamically to the presence of the solute and to its conformational changes. The coupling between the fluid and the protein allows reproducing accurately the reorganisation of the velocity field around the embedded protein, as illustrated in Supporting Information (SI) Figure S1, thus ensuring high physical realism.

6 ACS Paragon Plus Environment

Page 7 of 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Systems The unfolding process under shear was investigated for three systems: the peptide β-hairpin (fragment 41-56 at the C-terminal) from the GB1 protein (PDB code 1GB1), and two small proteins, the WW domain of 37 residues (PDB codes 1E0L) and a calcium binding domain constituted by 75 residues (PDB 1CLB). The secondary structures of the three proteins is represented in Figure 1. The three systems have been already used to test the performance of the LBMD technique based on the OPEP force field 25,29 . In this study, we used the version v.4 of the OPEP force field 37 .

Results Kinetics of unfolding We have initially investigated the characteristic time scale of the unfolding event at different shear rates. Unfolding was monitored by looking at the root mean square displacement (RMSD) of the protein backbone atoms with respect to their initial configuration (native state). For each unfolding event we recorded the first passage time (FPT) for the protein to unfold. FPT is defined as the time taken to observe a jump in RMSD larger than the threshold of 7 ˚ A. At each shear rate, 5 independent trajectories were generated for each protein in order to appreciate the degree of randomness of the process, and estimate the statistical error. The results are summarised in Figure 2. In the lower part of the figure, we plot a representative unfolding trace of the RMSD vs time for each system. The dividing cut-off to distinguish between folded and unfolded structures was set to 7 ˚ A, a value that captures well the systematic and abrupt jump in RMSD occurring during unfolding for all three proteins, see also SI Figure S2 where the probability distribution of RMSD for each system has been recorded at different temperatures so to distinguish the values associated to folded and unfolded states. In the upper part of

7 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 23

Figure 2 the mean FPT is reported as a function of the shear rate γ. ˙ For all proteins we observe two clear regimes. At highest shear rates, rapid unfolding occurs in a few nanoseconds time span and the kinetics weakly depends on γ. ˙ For the β-hairpin, the transition to the second regime where kinetics strongly depends on shear rate, corresponds to γ˙ ∼ 1011 s−1 , while for the two larger proteins, the WW and Ca-binding domains, the transition takes place at around 6 · 1010 s−1 . When decreasing γ˙ ∼ 3 · 1010 s−1 unfolding of β-hairpin and Ca-binding domain occurs in about one hundred nanoseconds. The WW domain still unfolds at very low value of γ˙ = 2.4 · 109 s−1 but again with a long mean FPT. At such low shear rate β-hairpin is stable and does not unfold during the entire simulation time of 200 ns. For this system we have also verified that the shear rate dependent kinetics is not affected by the choice of the RMSD cutoff (from 5 to 7 ˚ A) to distinguish the unfolding event, see Figure S3 in SI. Moreover, in order to check possible size effect, for the larger systems, WW and Ca-binding domains, we have performed extra calculations with a simulation box of doubled size. The unfolding kinetics is substantially the same, see Figure S4 in SI. Interestingly, the analysis of the mechanism of protein unfolding by force pulling has already shown the presence of two types of kinetic regimes depending on the magnitude of the applied pulling force 38–40 . The weak force-dependent regime corresponds to values of the force high enough to remove the unfolding barrier and making the transition dominated by the determistic dynamics of the system, while the force-dependent kinetic exponentially depends on the force magnitude with the unfolding rate-limited by the barrier, and thus is dominated by stochasticity 40 .

Structural weakness In order to probe the protein robustness, we have initially inspected all trajectories in order to identify the protein regions that are more prone to disruption by shear, that is, to identify the structural weak spots. A representative view is given in Figure 3 for all systems. Visually, the 8 ACS Paragon Plus Environment

Page 9 of 23

1

kMFPT (1/ns)

β-hairpin WW-domain Ca-binding 0.1

0.01

RMSD (Å)

0.1

RMSD (Å)

1 . -1 10 γ (s ) x10

10

15 10 5 0 0

1

2

3

4

5

6

7

8

9

10

1

2

3

4

5

6

7

8

9

10

1

2

3

4

5 6 time (ns)

7

8

9

10

15 10 5 0 0

RMSD (Å)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

15 10 5 0 0

Figure 2: Top Panel. Kinetics of the unfolding as a function of the shear rate γ˙ computed from the mean first passage time for unfolding. Error bars are estimated over five independent trajectories. Bottom Panel. Example of the time evolution of the RMSD for the three proteins. The unfolding time is obtained when RMSD becomes larger of a fixed threshold (7 ˚ A). unfolding of the β-hairpin peptide corresponds to the unzip process with progressive opening of the β-strands starting from the termini. For the other proteins, a clear weakness is individuated: in the WW domain, the three β strands forming the rigid core of the native structure are disassembled, possibly as a result of the elongational component of the shear flow. In particular the β3 strand at the C-terminal of the protein appears to be the weakest region. This strand disconnects from the rest of the protein and is stretched apart. Depending on the extension of the trajectory, all secondary structures get eventually disassembled. The resilience of the β1 -β2 motif agrees well with earlier investigations that demonstrated how the loop connecting β1 -β2 is the major stability factor for this fold 41 . For the Ca-binding domain, which native structure is represented by three connected α-helices, the unfolding mainly occurs at the N-terminal, with one of the helices (α1 ) disconnecting as structured element from the rest of the protein and 9 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 23

β-hairpin

WW domain β2 β1

β1 β3

β1 β2

β2

Ca domain α3

α3 α2 α1

α2

α1

α2

α1

α3 time

Figure 3: Characteristic structures of each protein along the unfolding path and extracted from a single simulation at γ˙ = 8·1010 . For the WW and the Ca binding domains, the labels indicate the secondary motifs in the protein fold. being stretched by the elongational forces. The initial separation of structured fragments has also been observed using the same type of simulations during the early step of unfolding of the A2 domain from the von Willibrand factor protein 26 . It was proposed by other authors that the same molecular mechanism for shear-induced unfolding applies to other large proteins 42,43 . The creation of protrusions is triggered by the periodic exposure to the elongational forces that causes the expansion and contraction of the globular shapes until a major instability triggers the disconnection of secondary structural elements.

Unfolding paths Sampling the unfolding path requires gathering meaningful statistics that is not always reachable by current computational means, therefore we focused on the simplest and smallest system, the β-hairpin, and performed a systematic investigation of the unfolding kinetics and its path-

10 ACS Paragon Plus Environment

Page 11 of 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

ways. We generated 50 independent initial conditions of the peptide in the absence of shear, these configurations were used as starting points for the runs at different shear rates. We first obtained a more complete vision of the unfolding kinetics by recording the distribution of mean FPTs for γ˙ ∈ [7 · 1010 − 1.3 · 1011 ]. The distributions of the unfolding time, as reported in SI Figure S2, features a long tail, and even at the highest shear rates, some unfolding events occur after almost a tens of nanoseconds. Since the thermal unfolding process of β-hairpin has been well characterised via enhanced sampling based on the OPEP force field 29 , we use here the same reaction coordinates that characterise the free energy landscape of the peptide. Namely, we computed the number of H-bonds formed between the strands, (HB) and the gyration radius of the hydrophobic core of the peptide (Rg ). The first collective variable is numerically evaluated by performing a smooth P 1−(d/d0 )6 counting of the HB pairs, HB = i,j 1−(d/d0 )12 , with the sum running over the distance d between the backbone amide hydrogens and carbonyl oxygens of residues i and j with j > i + 4 and setting the cut-off distance at d0 = 3˚ A. The second collective variable is taken by considering the heavy-atoms and the side-chains of the hydrophobic residues Trp3, Tyr5, Phe12 and Val14. We collected the traces of the unfolding trajectories in a time window of 600 ps centered at the unfolding time event [mean FPT-300ps, mean FPT+300ps], and projected them in the 2D space describing the free energy landscape of the peptide, as shown in Fig. 4. The unfolding events proceed along a clear main path (white dashed line), with the breaking of the HB spine anticipating the full opening of the peptide, and substantially reproducing, not surprising for this simple system, the path of the thermal unfolding previously studied 29 . To be noted that a handful of events proceeds along an alternative path (black line) where the system explores intermediate states featuring partial open structures.

11 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 23

U

I2 I1 F Figure 4: Characteristic unfolding path of the β-hairpin peptide monitored for 50 independent trajectories (each represented by a different coloured line), and projected in the 2D space defined by the number of backbone HB formed between the two strands, and the size of the hydrophobic core of the peptide (RG ).

The effect of drag force For many systems where shear flow triggers proteins conformational changes (e.g. catch bond systems) 3 , single molecule pulling experiments are routinely carried out to investigate the system response to external mechanical forces. In order to mimic the kind of conditions where unfolding is caused by a directional force acting along the end-to-end distance 44 , we performed steered MD simulations using a directional force of magnitude comprised between 50 pN and 500 pN . In our in silico essay on the β-hairpin, the first heavy atoms at one terminal of the peptide are constrained in space while the directional force (e.g. fx ) is applied at the opposite terminal end, see right panel in Figure 5. For each force magnitude, steered simulations were performed starting from 10 independent initial configurations. The unfolding kinetics was computed considering the extension of the end-to-end distance, and a cut-off value 6 ˚ A for discriminating between folded and unfolded configurations. The mean FPT for unfolding versus the magnitude of the pulling force is reported Figure 6. For very high force magnitudes (300 pN − 500 pN ) unfolding occurs in the sub-nanosecond timescale and the kinetics weakly depends on the magnitude 38–40 , similarly to what observed at the very high shear rates γ˙ ∼

12 ACS Paragon Plus Environment

Page 13 of 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Figure 5: Top. Molecular view of the perturbing forces acting on the β-hairpin under shear and under pulling experiments. Left panel: ramp of HBs creating the link between the two strands. Central panel: shear forces acting on the protein structure. Right panel: in silico pulling experiment, with a directional force acting on one atom of a HB pair (the other atom being constrained in space). 1011 s−1 . In the 50 pN − 200 pN regime, we recover the behaviour predicted by Bell’s model 44 that relates the rate of unfolding to the perturbing force, encoded by the expression:

−1 τunf =

kb T −β(∆G∗ −fx δx) e h

(2)

where the prefactor is given by the unidimensional transition state theory, (kB is the Boltzmann constant, T temperature and h the Plank constant), ∆G∗ is the unfolding free energy barrier along the separation reaction coordinate, and δx is a parameter that relates to the transition state along the reaction coordinate. By fitting our data, we obtain δx ' 1.5 ˚ A, a value that correlates with the fact that the unfolding is triggered by the separation of the HBs in the ramp. In addition, the fit provides ∆G∗ ' 8 kcal/mol, comparable to the energy of a few HBs. By comparing the unfolding kinetics for the peptide under shear and under pulling conditions, within our mechanical model for β-hairpin, we found that in the range γ˙ ' 109 − 1010 s−1 shear conditions are equivalent to a pulling force of 50 − 150 pN While the parallel with AFM-like pulling experiments is conceptually insightful, the effect of the shear flow has its own specificities. Under Couette flow, the peptide constantly tumbles due to the rotational component of the velocity field, so that its weak spots are cyclically exposed to the elongation forces, the second component of the shear field. In this way, the

13 ACS Paragon Plus Environment

The Journal of Physical Chemistry

10

kMFPT(1/ns)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 23

1 0.1 0.01 0.001 0

100

200 300 Force (pN)

400

500

Figure 6: Unfolding kinetics (1/τunf ) for the β-hairpin under the applied pulling force fx (pN ). The dashed line represents a fit according to Bell’s model. tension arising from the fluid drag affects simultaneously different regions of the molecule, and the magnitude of the shear forces on two portions of the molecule depends on their mutual distance. These aspects have been directly monitored in time during the unfolding events. According to a simple view, a static peptide like the β-hairpin would experience the maximal shear tension along the main unfolding reaction coordinate when the HBs ramp is well aligned with the main eigenvector of the shear field. However, caused by the tumbling peptide and the local reorganisation of the fluid, the magnitude of the shear forces along the unfolding reaction coordinate is highly modulated in time. At selected shear rates and initial configurations, we D ˜i for the HB donor (d) and acceptor have recorded the fluid component of the force FD i , fi = Γu

(a) backbone atoms (i) along the β strand. This contribution is then projected on the reaction coordinate defined as the HB separation distance (dHB ): ffrc = fdD · dHB − faD · dHB , see central panel in Figure 5. This has been done for all HB pairs along the folded β-hairpin. In the absence of shear, the average fluid perturbation < ffrc > for all HB pairs is practically zero. On the contrary, at high shear rate (γ˙ ∼ 1011 s−1 ), and before unfolding occurs, < ffrc >∼ 2 − 3 pN with instantaneous fluctuations as high as ∼ 40 pN . Moreover, the frequency of such high fluctuations increases with the magnitude of the shear rate, so that the higher the shear the more frequent and more intense is the impulsive force acting on the peptide. The frequency 14 ACS Paragon Plus Environment

20

dHB rc

10

ff cos(θ)

rc

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

dHB (Å) | ff (a.u.) | cos(θ)

Page 15 of 23

0

-10 0

1

2 time (ns)

3

4

Figure 7: Time evolution of the HB distance d corresponding to the first HB in the ramp at the C- and N termini (black line), the shear force projected on the HB ffrc (red), and the angular component cosθ = dHB · g where g is the shear gradient direction (green). The magnitude of the force ffrc has been scaled by an abritrary factor for the sake of clarity. of the shear forces correlates well with the characteristic time of the rotational component of the fluid that imposes the periodic alignment of the HB ramp with the velocity gradient with a characteristic time τR ∝ γ˙ −1 (e.g. see SI Figure S3). An overview of these aspects is presented in Figure 7 for one specific unfolding event (γ˙ ' 1011 s−1 ). It is interesting to notice that, at unfolding time, when the HB stretches so to reach the transition state, the tensile force increases because of this elongation. All this considered, in order to transpose the Bell-like model to shear perturbation, some ingredients must be plugged in. The external force is not constant in time but changes because of the rotational component of the velocity field, and secondly it increases along the reaction coordinate, fx (x, γ) ˙ = f0 (x, γ) ˙ · cos(γt), ˙ therefore on top of the randomness of the unfolding event, the shear flow imposes a time modulation of the kinetic barrier that depends on the rotational and the elongation dynamics of the molecule as a whole. While for a small peptide like β-hairpin the main effect of the tumbling motion is to change the exposure of the reaction coordinate, the HB ramp, to the shear stress, for a larger globular protein, tumbling induces cyclic deformations of the shape that can be followed as oscillation

15 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 23

of the protein gyration radius. In our case, shear-induced fluctuations are in the order of 20%, while for standard thermal fluctuation they are ∼ 10%. The oscillation ends with the critical separation, and disruption, of the weakest elements of the protein, leading to the formation of a protrusion that because of its separation from the body of the protein is more easily dragged by the shear, eventually leading to complete unfolding, as anticipated in section

42,43

.

Conclusions In this work we have explored the unfolding process triggered by shear flow on three model proteins of different size and fold. The simulations were based on the LBMD technique 25 , that couples a coarse-grained force field with the kinetic representation of the solvent, inclusive of hydrodynamic interactions. Namely, we exploited the coarse-grained OPEP force field 24 to characterise the proteins unfolding paths and timescales as a function of the shear rate. Our results unequivocally show that for the small proteins considered in this work shear-induced unfolding effectively takes place but only at high shear rates 5 . For computational feasibility, unfolding events were recorded at the timescale of tens of nanoseconds and occurred at shear rate of γ˙ > 109 s−1 . This corresponds to a shear strain for the unfolding time, W = γτ ˙ ≈ 102 . In principle, it is possible to extend our results to lower shear rates and longer timescales, so to target the critical value proposed by Jaspe and Hagen (γ˙ ∼ 107 s−1 ) for the unfolding of a small globular protein 5 . Such regime is today computationally expensive due to the requirement for sufficient statistics. In fact, assuming that several exposures to the elongation component are required before unfolding occurs, the characteristic time of the rotational field τR ' γ˙ −1 is on the order of 102 ns, while the unfolding event would occurs at the µs timescale. This, along with the focus on larger protein models, i.e. BSA 6 or even some units of the vWf 1,2 , would complete our understanding of shear-induced unfolding and will be the object of a subsequent investigation. Moreover, as experimental evidence suggests 10–12 , the conformational elongation due to shear

16 ACS Paragon Plus Environment

Page 17 of 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

forces could be critical for pathogen amyloid aggregation in physiological conditions. Our computational machinery lends itself to be a preferential tool to explore such effects. Thanks to the high resolution that OPEP reserves to the protein backbone, our approach allows to individuate the weakest regions of the proteins, and the simulations confirm the idea that the creation of a first protrusion triggers the complete unfolding of the molecule 43 . For the simple β-hairpin we demonstrated that unfolding under shear is similar to thermal unfolding and can also parallel atomic force microscopy experiments, where the termini of the peptide are pulled apart. The high shear rate (γ˙ ∼ 109 − 1010 s−1 ) needed to unfold the peptide in the nanosecond timescale (10 − 100 ns) corresponds to the in silico value for the directional mechanical force acting in pulling experiments in the order of 50 − 150 pN . However, we have shown that the tensile forces induced by shear, not only are cyclic and time dependent due the rotational component of the shear field, but also change along the reaction coordinate due to the different extension of the termini during the unfolding process.

Supporting Information Supporting data for the kinetics of unfolding (RMSD cut-off, size effects, large statistics, rotational dynamics).

Acknowledgement The research leading to these results has received funding from the ERC (FP7/2007-2013) Grant Agreement no.258748. Part of this work was performed using HPC resources from GENCI [CINES and TGCC] (Grant x201676818). We acknowledge the financial support by the ”Initiative d’Excellence” program from the French State (Grant ”DYNAMO”, ANR-11LABX-0011-01)”.

17 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 23

References (1) Springer, T. A. von Willebrand Factor, Jedi Knight of the Bloodstream. Blood 2014, 124, 1412–1425. (2) Fu, H.; Jiang, Y.; Yang, D.; Scheiflinger, F.; Wong, W. P.; Springer, T. A. Flow-Induced Elongation of von Willebrand factor Precedes Tension-Dependent Activation. Nature Comm. 2017, 8, 324. (3) Thomas, W. E.; Vogel, V.; Sokurenko, E. Biophysics of Catch Bonds. Ann. Rev. Biophys. 2008, 37, 399–416. (4) Kalas, V.; Pinkner, J. S.; Hannan, T. J.; Hibbing, M. E.; Dodson, K. W.; Holehouse, A. S.; Zhang, H.; Tolia, N. H.; Gross, M. L.; Pappu, R. V. et al. Evolutionary Fine-Tuning of Conformational Ensembles in FimH During Host-Pathogen Interactions. Science Advances 2017, 3, e1601944. (5) Jaspe, J.; Hagen, S. J. Do Protein Molecules Unfold in a Simple Shear Flow? Biophysical J. 2006, 91, 3415–3424. (6) Bekard, I. B.; Asimakis, P.; Teoh, C. L.; Ryan, T.; Howlett, G. J.; Bertolini, J.; Dunstan, D. E. Bovine Serum Albumin Unfolds in Couette Flow. Soft Matter 2012, 8, 385– 389. (7) Bekard, I. B.; Dunstan, D. E. Shear-Induced Deformation of Bovine Insulin in Couette Flow. J. Phys. Chem. B 2009, 113, 8453–8457. (8) Bekard, I. B.; Barnham, K. J.; White, L. R.; Dunstan, D. E. a-Helix Unfolding in Simple Shear Flow. Soft Matter 2011, 7, 203–210. (9) Lorna Ashton, . J. D.; Imomoh, E.; Balabani, S.; Blanch, E. W. Shear-Induced Unfolding of Lysozyme Monitored In Situ. Biophysical J. 2009, 96, 4231–4236. 18 ACS Paragon Plus Environment

Page 19 of 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

(10) Hamilton-Brown, P.; Bekard, I.; Ducker, W. A.; Dunstan, D. E. How Does Shear Affect A Fibrillogenesis? J. Phys. Chem. B 2008, 112, 16249–16252. (11) Trumbore, C. N. Shear-Induced Amyloid Formation in the Brain: I. Potential Vascular and Parenchymal Processes Shear-Induced Amyloid Formation in the Brain: I. Potential Vascular and Parenchymal Processes Shear-Induced Amyloid Formation in the Brain: I. Potential Vascular and Parenchymal Processes. J. Alzheimers Dis. 2016, 54, 457–470. (12) Byington, M. C.; Safari, M. S.; Conrad, J. C.; Vekilov, P. G. Protein Conformational Flexibility Enables the Formation of Dense Liquid Clusters: Tests Using Solution Shear. J. Phys. Chem. Lett. 2016, 7, 2339–2345. (13) Br¨ uckl, L.; Schr¨oder, T.; Scheler, S.; Hahn, R.; Sonderegger, C. The Effect of Shear on the Structural Conformation of rhGH and IgG1 in Free Solution. J. Pharm. J. 2016, 105, 1810–1818. (14) Schneider, S. W.; Nuschele, S.; Wixforth, A.; Gorzelanny, C.; Alexander-Katz, A.; Netz, R. R.; Schneider, M. F. Shear-Induced Unfolding Triggers Adhesion of von Willebrand Factor Fibers. Proc. Natl. Acad. Sci. USA 2007, 104, 7899–7903. (15) Sing, C. E.; Alexander-Katz, A. Elongational Flow Induces the Unfolding of von Willebrand Factor at Physiological Flow Rates. Biophysical J. 2010, 98, 35–37. (16) Szymczak, P.; Cieplak, M. Stretching of Proteins in a Uniform Flow. J. Chem. Phys. 2006, 125, 164903. (17) Szymczak, P.; Cieplak, M. Protein in Shear Flow. J. Chem. Phys. 2007, 127, 155106. (18) Alexander S. Lemak, J. R. L.; Chen, J. Z. Y. Unfolding Proteins in an External Field: Can We Always Observe the Intermediate States? Phys. Rev. E 2003, 67, 031910.

19 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 23

(19) Taketomi, H.; Ueda, Y.; G¯o, N. Studies on Protein Folding, Unfolding and Fluctuations by Computer Simulation. I. The Effect of Specific Amino Acid Sequence Represented by Specific Inter-Unit Interactions. Int. J. Pept. Prot. Res. 1975, 7, 445–459. (20) Kaya, H.; Chan, H. S. Towards a Consistent Modelling of Protein Thermodynamic and Kinetic Coperativity : How Applicable is the Transition State Picture to Folding and Unfolding? J. Mol. Biol. 2002, 315, 899–909. (21) P. Das, S. M.; Clementi, C. Balancing Energy and Entropy: a Minimalist Model for the Characterisation of Protein Folding Landscapes. Proc. Natl. Acad. Sci. USA 2005, 102, 10141–10146. (22) Chen, Z.; Lou, J.; Zhu, C.; Schulten, K. Flow-Induced Structural Transition in the bSwitch Region of Glycoprotein Ib. Biophys. J. 2008, 95, 1303–1313. (23) Lou, J.; Zhu, C. Flow Induces Loop-to-b-Hairpin Transition on the b-Switch of Platelet Glycoprotein Iba. Proc. Natl. Acad. Sci. USA 2008, 105, 13847–13852. (24) Sterpone, F.; Melchionna, S.; Tuffery, P.; Pasquali, S.; Mousseau, N.; Cragnolini, T.; Chebaro, Y.; St-Pierre, J.-F.; Kalimeri, M.; Barducci, A. et al. The OPEP Protein Model: from Single Molecules, Amyloid Formation, Crowding and Hydrodynamics to DNA/RNA Systems. Chem. Soc. Rev. 2014, 43, 4871–4893. (25) Sterpone, F.; Derreumaux, P.; Melchionna, S. Protein Simulations in Fluids: Coupling the OPEP Coarse-Grained Force Field with Hydrodynamics. J. Chem. Theory Comput. 2015, 11, 1843–1853. (26) Chiricotto, M.; Sterpone, F.; Derreumaux, P.; Melchionna, S. Multiscale Simulation of Molecular Processes in Cellular Environments. Philos. Trans. A 2016, 374, 20160225.

20 ACS Paragon Plus Environment

Page 21 of 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

(27) Sterpone, F.; Nguyen, P.; ing

the

Doutreligne, S.;

Tran, T. T.;

Melchionna, S.;

Baaden, M.;

Derreumaux, P. Multi-Scale Simulations of Biological Systems Us-

OPEP

Coarse-Grained

Model.

Bioch.

Biophys.

Res.

Comm.

2017,

https://doi.org/10.1016/j.bbrc.2017.08.165 . (28) Barducci, A.; Bonomi, M.; Derreumaux, P. Assessing the Quality of the OPEP CoarseGrained Force Field. J. Chem. Theory Comput. 2011, 7, 1928–1934. (29) Sterpone, F.; Nguyen, P.; Kalimeri, M.; Derreumaux, P. Importance of the Ion-Pair Interactions in the OPEP Coarse-Grained Force Field: Parametrization and Validation. J. Chem. Theory. Comput. 2013, 9, 4574–4584. (30) Kalimeri, M.; Derreumaux, P.; Sterpone, F. Are Coarse-Grained Models Apt to Detect Protein Thermal Stability? The Case of OPEP Force Field. J. Non Crys. Sol. 2015, 407, 494–501. (31) Ahlrichs, P.; D¨ unweg, B. Lattice-Boltzmann Simulation of Polymer-Solvent Systems. Int. J. Mod. Phy. C 1998, 9, 1429–1438. (32) Ahlrichs, P.; D¨ unweg, B. Simulation of a Single Polymer Chain in Solution by Combining Lattice Boltzmann and Molecular Dynamics. J. Chem. Phys. 1999, 111, 8225–8239. (33) Limbach, H.; Arnold, A.; Mann, B.; Holm, C. ESPResSo—an Extensible Simulation Package for Research on Soft Matter Systems. Comput. Phys. Comm. 2006, 174, 704–727. (34) Succi, S. The Lattice Boltmzann Equation for Fluid Dynamics and Beyond ; Clarendon Press, 2001. (35) Benzi, R.; Succi, S.; Vergassola, M. The Lattice Boltzmann Equation: Theory and Applications. Phys. Rep. 1992, 222, 145–197.

21 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 23

(36) Chiricotto, M.; Melchionna, S.; Derreumaux, P.; Sterpone, F. Hydrodynamic Effects on b-Amyloid (16-22) Peptide Aggregation. J. Chem. Phys. 2016, 145, 035102. (37) Chebaro, Y.; Pasquali, S.; Derreumaux, P. The Coarse-Grained OPEP Force Field for Non-Amyloid and Amyloid Proteins. J. Phys. Chem. B 2012, 116, 8741–8752. (38) Szymczak, P.; Cieplak, M. Stretching of Proteins in a Force-Clamp. J. Phys.: Condens. Matter 2006, 18, L21. (39) West, D. K.; Brockwell, D. J.; Olmsted, P. D.; Radford, S. E.; Paci, E. Mechanical Resistance of Proteins Explained Using Simple Molecular Models. Biophys. J. 2006, 90, 287–297. (40) Sun, L.; Noel, J. K.; Sulkowska, J. I.; Levine, H.; Onuchic, J. N. Connecting Thermal and Mechanical Protein (Un)folding Landscapes. Biophys. J. 2014, 107, 2950–2961. (41) Nguyen, H.; J¨ager, M.; Moretto, A.; Gruebele, M.; Kelly, J. W. Tuning the Free-Energy Landscape of a WW Domain by Temperature, Mutation, and Truncation. Proc. Natl. Acad. Sci. USA 2003, 100, 3948–3953. (42) Alexander-Katz, A.; Schneider, M. F.; Schneider, S. W.; Wixforth, A.; Netz, Shear-FlowInduced Unfolding of Polymeric Globules. Phys. Rev. Lett. 2006, 97, 138101. (43) Alexander-Katz, A.; Netz, R. R. Dynamics and Instabilities of Collapsed Polymers in Shear Flow. Macromolecules 2008, 41, 3363–3374. (44) Evans, E.; Ritchie, K. Dynamic Strength of Molecular Adhesion Bonds. Biophys. J. 1997, 72, 1541–1555.

22 ACS Paragon Plus Environment

Page 23 of 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

TOC 71x47mm (72 x 72 DPI)

ACS Paragon Plus Environment