How Fast Is Too Fast in Force-Probe Molecular Dynamics Simulations

Subscriber access provided by BUFFALO STATE

Article

How Fast Is Too Fast in Force-Probe Molecular Dynamics Simulations? Steven Sheridan, Frauke Gräter, and Csaba Daday J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.9b01251 • Publication Date (Web): 10 Apr 2019 Downloaded from http://pubs.acs.org on April 14, 2019

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 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 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.

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 12 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

How Fast Is Too Fast in Force-Probe Molecular Dynamics Simulations? Steven Sheridan1 , Frauke Gräter1,2,* , and Csaba Daday1,** 1

Heidelberg Institute for Theoretical Studies, Schloß-Wolfsbrunnenweg 35, 69118 Heidelberg, Germany 2 Interdisciplinary Center for Scientific Computing (IWR), Heidelberg University, Mathematikon, INF 205, 69120 Heidelberg, Germany * [email protected] ** [email protected] March 26, 2019

Abstract

5-100 m/s, which is about two orders of magnitude slower than the expected speed of sound. Our results can help in identifying and removing this bias While Molecular Dynamics (MD) simulations are from future simulations. routinely used to interpret atomic force microscopy (AFM) experiments of protein unfolding, computational cost in MD simulations still mostly imposes Introduction a large difference in loading rates and time scales in 1 this comparison. Loading rate dependencies of unfolding forces and mechanisms have been studied in Single molecule force spectroscopy is an established depth in experiments, simulations and theory. One technique to study biomolecules under pulling potential additional implication of the larger MD forces. Molecular Dynamics simulations are roupulling velocity that remains to be assessed is that tinely used to interpret these force spectroscopy regions of the proteins that are close to the point of experiments, but, due to computational restricforce application will be under force earlier or under tions, typically use pulling velocities on the order more force than more shielded regions, resulting in of 10−3 m/s to 1 m/s, as compared to those used a bias of the protein unfolding sequence which is in force spectroscopy experiments, approximately likely marginal at the slower AFM velocities. We 10−5 m/s to 10−8 m/s. This difference of several here, for the first time, quantify the parameters of orders of magnitude also holds for the applied loadthis bias using a model system of four tandem spec- ing rates, the product of the pulling velocity and trin repeats (SRs) linked with long, flexible poly- the system’s spring constant. Recently, high-speed glycine linkers. We subject the system to seven AFM has reached MD pulling velocities and loading different pulling velocities ranging from 0.01 to 10 rates, resulting in a good agreement between the m/s and find that for the fastest velocities, down experiments and simulations in terms of unfolding to 1 m/s, the outer domains preferentially unfold; forces. 1,2 One alternative method that can attempt in fact, at 10 m/s, this happened in 100 cases out to bridge this gap is through Boxed Molecular Dyof 100. Based on these data, and also through an- namics (BXD) 3 although at the cost of losing exact alyzing the amount of partial unfolding in the be- time scale information and additional approximaginning of the simulations, we show that the bias tions. is equivalent to an effective signal propagation of The effect of the loading rate on the force re1

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

sponse of biomolecules has been investigated at various levels of detail, in theory and experiments. 4–6 The order of unfolding can vary with the loading rate as the molecule can cross a multitude barriers on the multidimensional energy landscape of the complex, which are differently affected by the external loading. 7 An additional complexity is given by the finite time that is required for an external force to propagate through the molecule, a question not taken into account by the models described above. Speeds of force propagation or, more generally, signal propagation have been obtained using a diverse set of techniques and systems. Low-frequency modes in vibrational spectra revealed a speed of sound of 2300 m/s in myoglobin. 8 This order of magnitude is in line with the 300-500 fs time scale for signal propagation in myoglobin (corresponding to 3500-5800 m/s speed of sound assuming a 3.5 nm protein diameter) observed in time-resolved serial femtosecond crystallography 9 and combined smallangle X-ray scattering experiments and MD simulations. 10 We have observed in two very different systems, an unfolded peptide chain and a lipid bilayer, velocities of force propagation of 5000 m/s and 1400 m/s, respectively, using MD simulations. 11,12 Taken together, the speed of sound of proteins is in line with the values expected for soft matter systems and liquids in general (e.g. the speed of sound of water is 1500 m/s) 13 ). The finite time of an external force to propagate through a protein or even a multi-domain protein of up to 10 nm diameter should in this light not be larger than a few picoseconds. The known speed of sound in proteins would thus argue against a role of the force propagation speed in the unfolding sequence and mechanism of proteins, as on the time scale of unfolding, the short delay in loading for some inner domains or subdomains relative to those domains close to the force application points would be marginal. However, the mode of force propagation and relaxation relevant for force-induced protein unfolding, including the role of friction known to affect rupture forces at high MD pulling speeds, 14 is unknown. Molecular Dynamics simulations often yield force-induced unfolding pathways starting with the unfolding from the termini, also for multi-domain proteins with domains of similar mechanical stability. 15 It is an open question how the unfolding probability of a

domain or subdomain within the protein under investigation is affected by its proximity to the termini of the construct, where the external perturbation is applied. We sought to determine whether proximity of a (sub)domain to the point of force application biases the likelihood of the domain unfolding, and how any such effect depends on the loading rate. To this end, we performed force-probe molecular dynamics simulations of a model system, based on spectrin. Spectrin repeats are alpha-helical protein domains the mechanical unfolding of which has been studied in depth in both single molecule force spectroscopy experiments and MD simulations. 16,17 Our construct comprises four spectrin domains with identical sequence, which allows direct quantification of the bias from the deviation of equal probabilities. Including large numbers of identical domains is also commonplace in force spectroscopy experiments. This ensures that each single approach-retract cycle captures several unfolding events, thereby increasing the throughput of the experiments. We separated the domains by a long (more than 3 nm when fully stretched) interdomain linker to exclude cooperativity between the domains. We assessed a wide range of velocities typical for force-probe molecular dynamics, spanning three orders of magnitude, from 0.01 m/s through 10 m/s. We found that at pulling rates 1 m/s and higher, the order of unfolding was significantly biased toward the outer domains. We infer a speed of propagation of between 5-100 m/s through two independent analyses, the loading rate dependence of the bias and the time delay between the stretching of different repeats. This estimate is around 100 times below previously reported speeds of tension propagation through proteins. We found external friction and intramolecular rearrangements to both substantially contribute to the unexpectedly slow apparent force propagation. Our results can guide future MD simulations and can help to put them in context of single-molecule force spectroscopy experiments. More generally, they can help in understanding the force response of proteins.

2

Methods

Simulations were conducted using GROMACS 2016, 18 with the Amber 99SB* force field 2


Page 2 of 12

Page 3 of 12 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


with ILDN modifications, 19,20 the TIP3P water model, 21 and Joung ions. 22 We used LINCS holonomic constraints 23 for all bonds, and used virtual sites 24 for all hydrogen atoms. Thanks to these constraints, we used 5fs time steps. This choice is validated against 2fs femtosecond time steps in the Supplementary Material. We employ a Bussi-Parrinello velocity rescaling thermostat 25 at 300 K and a time constant of 0.1 ps for both the protein and the non-protein atoms, and an isotropic Parrinello-Rahman barostat 26 with a time constant of 1 ps set to 1 bar and a compressibility of 4.5·10−5 bar−1 . Long-range electrostatic interactions were treated through a fourth-order Particle Mesh Ewald (PME) 27 method with a cutoff radius of 1 nm and a 0.16 nm grid spacing. A Verlet neighbor list was updated with a frequency tuned automatically by GROMACS but no more often than every 6 steps (30 fs). We created a construct of four repeats of the chicken brain alpha spectrin repeat (SR) 17, from the PDB structure 1CUN. 28 To acquire starting structures, we first performed 100 ns of equilibrium simulation on SR17 and extracted 10 equidistant frames in time starting at 10 ns. Then, we constructed a tandem structure of four identical copies of these frames linked with three linkers of 10 glycines each. This tandem structure was aligned with the x-axis of a box with dimensions 64×6.5×6.5 nm3 , solvated in water and Na, Cl ions at 0.1 M. The final system had approximately 270,000 atoms. For force-probe simulations, we applied a harmonic (umbrella) potential onto both termini of the construct, with a spring constant of 830 pN/nm. Pulling speeds were 10, 10/3, 1, 1/3, 1/10, 1/30, and 1/100 m/s. We did 100 simulations of 10, 3.3, 1, and 0.33 m/s; and 10 simulations each of 0.1, 0.033, and 0.01 m/s. We performed a total of 160 5 ns-long constant-force pulling simulations across a range of forces, from 21 to 332 pN, at steps of 21 pN, with 10 replicates at each force level. The results of these constant-force pulling simulations are discussed in the Supplementary Material. Furthermore, we repeated the equilibrium simulations and performed 10 pulling simulations at 1 m/s with a time step of 2 fs for validation. We obtained rupture forces from the force profiles as in: 29 we took the average of the two forces at

the termini and applied a Gaussian smoothing with a width consistent with an extension of 0.1 nm. We then took the highest value of this smoothed force out of the first 15 nm of total extension, which corresponds to the main unfolding event. For the determination of the loading rate, we took the highest difference between any point of the aforementioned smoothed force profile that were within the first 10 nm of total extension and the force at 0.5 nm. The smoothing is necessary as force profiles have high noise levels and a moving window is necessary as the approximately linear increase in the force profiles does not begin in the same point in all simulations. The following measures were collected for each repeat and each glycine linker at 1ps time resolution: length, solvent accessible surface area (SASA), and secondary structure (using DSSP 30 ). For the primary analysis, we used the end-to-end distance of individual spectrin repeat units.

3

Results

Force Propagation from Unfolding Bias We subjected the protein construct, four spectrin domains separated by short glycine linkers, to a total of 430 force-probe Molecular Dynamics simulations at four different speeds between 0.01m/s and 10m/s. In all of these simulations, after an initial extension from the straightening of the glycinelinkers (Figure 1), we observed the unfolding of individual spectrin repeats as distinct events: once a repeat began unfolding, it continued, and in most cases it finished before the next began. At 10 m/s and 3.3 m/s, however, domains occasionally unfolded simultaneously. To quantify whether some domains preferentially unfold relative to others, we monitored the average end-to-end distances of the constructs and the individual spectrin repeats therein. We present the preferential unfolding of the outer spectrin repeats in Table 1 and Figure 2. At high speeds (1 m/s and above), there was a strong, statistically significant preference for the outer spectrin repeats to unfold first. The bias was more evident with increasing pulling speeds, with 10m/s pulling velocities yielding exclusively an outer do3



Page 4 of 12

Relaxed +6 nm Stretched %inner

%outer Outer SR breaks

+18 nm

Inner SR breaks

Figure 1: The two different unfolding pathways of the system: after a flexible stretching and partial unfolding of the spectrin repeats, amounting to 6 nm of the unfolding, either one of the outer domains (blue) or the inner domains (green) unfolds. All three poly-glycine linkers are colored gray. speed (m/s) 10.00 3.333 1.000 0.333 0.100 0.033 0.010

# inner 0 14 34 48 4 4 6

# outer 100 86 66 52 6 6 4

a.

p-value 2·10 −30 8·10 −14 0.002 0.8 0.8 0.8 0.8

10 m/s

c.

Table 1: The dependence on the pulling speed of the preference for outer domains rather than inner ones to unfold. The p-value is computed from a binomial test.

b.

1 m/s

d.

Figure 2: The outcomes of the 100 simulations at 10 m/s (a) and 1 m/s (b); 10 random simulations are highlighted in each. The ratio of simulations that showed an outer domain (blue in a and b) fits on a logistic function of the loading rate (c). The apparent elastic constant of the system (kprot ) strongly depends on the pulling speed (d).

main spectrin unfolding as the first event. A speed of 0.33 m/s resulted in an insignificant difference between the number of inner and outer unfoldings among the 100 simulations, suggesting that only a very minor, if any, bias is present at this velocity in case of our test system. Interestingly, among the outer repeats, we also observed a preference for repeat 4 over repeat 1 at higher velocities (71-29 at 10 m/s), which was again lost below 1 m/s. This is understandable since the external perturbation is acting on a different helix of the repeat, underlining the anisotropic response of proteins with respect to the force application points. The dependence of the preference on the loading rates fit well on a logistic function. To first approximation, this can be understood using the Bell model: 4


Page 5 of 12 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


k1,4 k1,4 + k2,3

= =

transition state of 0.11 nm. The fact that the two different models led to opposing corrections to the Bell model’s value, and that both corrections alone yielded a nearly identical fit (see Fig. S3), leads us to conclude that in this case, the distance to the transition state obtained from the Bell model is an acceptable approximate value. We note that we determined the loading rates used in Fig. 2c and Fig. S3 directly from the force profiles, and therefore also were able to estimate the stiffness kprot of the construct. The tandem spectrin repeats including the glycine linkers feature an overall stiffness in the range of 15-75 pN/nm (Fig. 2d). Interestingly, this stiffness value was found to increase at higher loading rates, varying by a factor of nearly five, in line with a time-dependent, i.e. viscoelastic response of the system. The pulling velocity-stiffness relationship has been described to be logarithmic before, 31 although further analysis showed the elastic constant to tend to a constant value in the limit of low pulling velocities. 32 Our data have a slight convexity in the log-linear plot, but we will not attempt to analyze this further as it is not the main focus of our work.

e−∆G1,4 /kB T (1) e−∆G1,4 /kB T + e−∆G2,3 /kB T 1 1 = 1 + e∆∆g 1 + e−∆F ∆x/kB T

The advantage of the Bell model is that in this case, this bias (i.e., the ratio of the unfolding rates) is constant in time, meaning that we need not define a single moment of time for the rupture. Assuming that there is a time delay ∆t between outer and inner domains, we obtain the difference in barrier heights as: ∆∆g =

F˙ · ∆t · ∆x ∆F · ∆x = , kB T kB T

(2)

which in turn will mean that the loading ratedependence of the preferential unfolding is a logistic function: k1,4 1 1 . = = k1,4 + k2,3 1 + e∆∆g 1 + e−aF˙

(3)

In the above, ∆x is the distance to the transition state which we obtained from a Bell fit of the rupture forces (Figure S3). It here refers to an effective time delay that explains the bias in unfolding probabilities. A microscopic interpretation will follow further below. From the force profiles, we found the distance to the transition state, ∆x, to be 0.08 nm, and therefore the time delay between the forces experienced by the outer and the inner spectrin repeats, ∆t, to be 550 ps. From this ∆t and the equilibrium inter-domain distance of 7.6 nm, we could infer an apparent velocity of force propagation of 14 m/s. As expected from the several orders of magnitude in pulling speeds, the fit of the rupture force with the loading rate was unsatisfactory, as the rupture forces observed at high pulling speeds werere systematically higher than the ones expected from a simple extrapolation of lower pulling speeds as in the Bell model. The Bell model can be extended in two different ways. On the one hand, using the Dudko-Hummer-Szabo model, 4 accounting for the changing position of the transition state, we get a decreased distance to the transition state of 0.04 or 0.06 nm (for ν=1/2 or 2/3, respectively). On the other hand, by including a linear term for friction, 14 we obtained an increased distance to the

Force Propagation from Partial Unfoldings To further explore this speed of force propagation, we also examined the time series of subdomain extension. The first approximately 20% of each pulling simulation featured approximately linear force-extension curves. Each spectrin repeat and glycine linker contributed to this overall extension to a different extent. We monitored the individual extensions of each SR and glycine linker during the initial straightening of the protein construct (Fig. 3a, top). We then inferred the time delay ∆tSR by superimposing the time series of the extensions of the inner two spectrin repeats to those of the outer two spectrin repeats, and maximizing their overlap. Analogously, as a second independent measurement, we repeated this analysis for the extension of the central glycine linker and its time delay relative to the two outer glycine linkers, resulting in ∆tGLY . Fig. 3 shows the speed of force propagation obtained from the time delay in glycine and spectrin repeat extensions. Based on this partial unfolding analysis, we infer speeds of propagation between 8-84 m/s across the different pulling ve5



reference time (ps) speed (m/s) Yu&Leitner 8 2300 Barends et al 9 0.5 35001 Brinkmann&Hub 10 0.3 58001 11 Young et al 5000 Aponte et al 12 14002 Bilaniuk&Wong 13 15003 4 this work 550 14 this work5 90-1000 8-84 1: Assuming a diameter of myoglobin of 3.5 nm. 2: Based on simulations on lipid bilayers. 3: In pure water at 300 K. 4: Based on a constant time delay explaining unfolding preferences. 5: Based on observed partial unfoldings.

We estimate the friction coefficient of a single spectrin repeat in water in additional MD simulations. We pull the protein through water at a constant velocity ranging from 0.1-4 m/s while keeping the vector connecting its termini parallel with the pulling direction with another restraint force, thereby mimicking the alignment of the domain along the pulling direction present within the larger construct. We observe the friction force to depend linearly on the pulling speed with a slope of 35 pN/(m/s) (Fig. S2). The friction coefficient is within 10% of this value in the case of only one restraint force. We next compare this friction acting on an individual spectrin domain in water to the apparent friction that is required to explain the bias in domain unfolding. To this end, we can re-formulate the sigmoidal dependence of the preferential unfolding in Equation 2 by substituting F˙ · ∆t by γ · v1,4 . We measure the average velocity of the center of mass of the outer domains and find them to be about ±40% of the pulling velocity (meaning that about ±10% is stretching of the spectrin repeats and/or the external springs). Using these velocity values, we can estimate the frictional coefficient needed to explain the preferential unfolding as γ =75 pN/(m/s), assuming a distance to the transition state of 0.08 nm. This means that friction with water explains a sizable portion, i.e., about 50%, but not all of the observed preferential unfolding. The remaining bias is most likely due to internal friction.

Table 2: Speeds of sound in proteins and other soft matter obtained previously using experimental or simulation techniques in comparison to the speeds of force propagation obtained here.

locities, with both measurements (glycine and SR extensions) yielding comparable speeds within the error. There is no apparent trend with pulling velocity. This range of speeds is well in line with the previous, independently determined 14 m/s. The spread of these values among individual trajectories is low at high pulling speeds, and high at lower pulling speeds. However, across the whole spread, the propagation speeds are substantially lower than the speed of sound in proteins and other soft matter, which is approximately 1000-5000 m/s (Table 2, Figure 3b). We obtain similar time delays for 2 fs time steps (370-630 ps) and constant-force pulling Intramolecular Friction Contributes simulations (130-410 ps) as well, as shown in Figure to the Unfolding Bias S2. Since our spectrin repeats are helical domains, internal friction can be reasonably interpreted as Friction with Water Only Partially inter-helix rearrangement. We assess the extent of internal conformational rearrangements by analyzExplains the Unfolding Bias ing if the overall increase in spectrin repeat extenOne explanation for the decreased force acting on sions under force can be solely explained by helix the inner as compared to the outer domains is fric- unwinding or if additional conformational changes tion with water. The outer domains transmit only add to the overall extension. To this end, we anala part of the external perturbation force, since yse the 100 1 m/s trajectories in terms of their endfriction with water slows down their motion and to-end distance and helical residues (Fig. 4a) We thereby the force transmitted to the inner domains observe a large spread in terms of the number of heis decreased. At higher pulling velocities, this can lical residues in the repeat, Nα , and the correspondbe a considerable effect slowing down protein un- ing absolute end-to-end distance, ∆x, suggesting folding. 2,14 that helix unwinding does not fully explain spec6


Page 6 of 12

Page 7 of 12 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


a.

Δt SR

a.

b.

b.

Δt Gly

Figure 4: A heatmap of extension and helicity of the 100 simulations performed at 1 m/s (a) shows a broad range of possible extensions. All four spectrin repeats are pooled together in this case. From the 400 different time series, we obtain a wide variety of helix/extension slopes, as can be seen in the histograms (b). Highlighted is the slope of -6.7 helical residues/nm, corresponding to an ideal helix unwinding into a worm-like chain.

Figure 3: (a) The speed of force propagations determined from partial unfoldings (symbols labeled with the pulling speeds in m/s) and from the logistic fit and the Bell model (labeled by “Bell”). (b) The values found in our study compared to the speed of sound reported in literature (orange symbols). The error bars show the speed determined by the mean time delay +/- one standard error on the mean. For the Bell model, the uncertainty stems bution to the overall ‘friction’ causing the observed from the distance to the transition state, 0.04-0.08 bias in domain unfolding. nm.

Discussion and Conclusions

trin repeat extension during unfolding. For each trajectory, we obtained the average number of lost helical residues per nanometer extension from the slope, ∆Nα /∆x, of a linear fit to this data. We find a wide range of outcomes with a mean slope of 3.9 residues/nm and an inter-quartile range of 2.6-5.2 residues/nm (Fig. 4b). For the unfolding of an ideal helix, one instead expects 6.7 residues/nm, or 0.15 nm/residue, resulting from the difference between the length per aminoacid of 0.3 nm in an unfolded peptide chain (at 100 pN, based on a worm-like chain behavior), and the rise in helix length of 0.15 nm per amino acid. Thus, the significantly lower number of residues we find to yield a nanometer of extension, 3.9 instead of 6.7 residues/nm means that on average, about 60% of the “flexible response” of spectrin repeats is due to helical unwinding and 40% is due to helical rearrangement. We note that this 40% is an average, and the balance between unwinding and rearrangement in response to the external force varies widely (Fig. 4b). We conclude that intradomain conformational changes contribute a signficant portion to the force-induced extension and are a likely contri-

In this work, we have shown on a model tandem protein that regions of proteins that are closer to external mechanical force are more likely to unfold than more shielded ones. Importantly, this is not due to cooperativity nor due to differences in intrinsic stability of the domains, as both of these factors were excluded in the setup of the system. While this is an expected property of these systems, the exact parameters of this “proximity bias” have not been quantified before. While these parameters will undoubtedly vary depending on the system, our work shows that the effective force propagation speed is about two orders of magnitude lower than the speed of sound, which is the most straightforward reference for comparing the propagation of mechanical signals. The surprisingly slow speed of propagation means that more care is needed than previously assumed when interpreting force propagation and preferential unfoldings. It should not be excluded that due to the high loading rates in simulations, a domain preferentially unfolds not because it is necessarily less stable than other domains, but because 7



it is closer to the point of force application. Our results suggest that a pulling speed below 1 m/s is advisable to avoid such artifacts, which again should be understood as an order of magnitude estimation which is expected to vary with the system under investigation. For example, a change in unfolding order was previously not observed 15 despite changing the oxidation state of domains when pulled at a velocity of 2.5 m/s, a result that is more understandable in light of our findings. While the increase in computing capacities and code efficiencies of MD have allowed a significant reduction of the MD pulling speeds, many studies still remain close to the limit of 1m/s we determined here. Also, friction with water and internal friction, among other determinants, scale with the size of the system. Hence, larger proteins or protein assemblies will already face an unfolding bias among domains or subdomains at lower velocities, which are moreover more difficult to reach for these systems. Simulating at multiple pulling speeds can help determine whether the speed of pulling is biasing results. Alternatively, if the protein can be approximated as a tandem structure, extrapolating preferential unfoldings from high pulling speeds to low ones can be feasible. In our case, the repeats were equally stable by construction. Instead, in realistic cases comprising domains different in structure and stability, the “true” preference between inner and outer regions can still be recovered from a logistic regression similar to the one presented in this work, even if the bias due to a high pulling speed was not fully eliminated in the systems. We find intramolecular rearrangements to be involved in extended spectrin under tension and to likely be a major component of the viscous response and slowed-down force propagation. While helix unwinding has been the major focus of previous studies of spectrin’s force response 17 our study emphasizes the role of inter-helical packing and orientation. In particular, roughly half of the SR extension is not explained by helix unwinding but other conformational changes, underlining the need for extensive sampling even in the case of this ostensibly simple topology. It remains to be investigated how other protein systems, from stiff immunoglobulin domains to soft ankyrin repeats, might differ in their intramolecular friction and force propagation speeds. Experimental loading rates and most physiologi-

cal conditions are well below the range of velocities considered here, and would be expected, from our results here, to be free of bias toward outer regions at length scales comparable to individual proteins. Importantly, however, our work shows that the relevant timescales of force relaxation are much longer than would be expected from a simple mechanical picture, a finding likely important for force propagation at mesoscopic length scales.

Supporting Information Available Water-spectrin drag force raw data, force-extension curves at constant force, loading rate-rupture force dependence, and the validation of the speed of propagation of the force signal.

Acknowledgments C. D. and F. G. acknowledge funding from the Klaus Tschira Foundation. The authors acknowledge support by the state of Baden-W¨ urttemberg through bwHPC and the German Research Foundation (DFG) through grant INST 35/1134-1 FUGG. C. D. thanks Helmut Grubm¨ uller and Carleen Kluger for helpful discussions.

References [1] Rico, F.; Gonzalez, L.; Casuso, I.; PuigVidal, M.; Scheuring, S. High-Speed Force Spectroscopy Unfolds Titin at the Velocity of Molecular Dynamics Simulations. Science 2013, 342, 741–743. [2] Rico, F.; Russek, A.; Gonzalez, L.; Grubm¨ uller, H.; Scheuring, S. Heterogeneous and Rate-Dependent Streptavidin-Biotin Unbinding Revealed by High-Speed Force Spectroscopy and Atomistic Simulations. Proc. Natl. Acad. Sci. U.S.A. (published ahead of print) [3] Booth, J. J.; Shalashilin, D. V. Fully Atomistic Simulations of Protein Unfolding in Low Speed Atomic Force Microscope and Force Clamp 8


Page 8 of 12

Page 9 of 12 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


[4]

[5]

[6]

[7]

[8]

Experiments with the Help of Boxed Molec- [15] Lee, E.; Hsin, J.; von Castelmur, E.; ular Dynamics. J. Phys. Chem. B 2016, 120 . Mayans, O.; Schulten, K. Tertiary and Secondary Structure Elasticity of a Six-Ig Titin Dudko, O. K.; Hummer, G.; Szabo, A. IntrinChain. Biophys. J. 2010, 98 . sic Rates and Activation Free Energies from Single-Molecule Pulling Experiments. Phys. [16] Randles, L. G.; Rounsevell, R. W.; Clarke, J. Rev. Lett. 2006, 96, 108101. Spectrin Domains Lose Cooperativity in Forced Unfolding. Biophys. J. 2007, 92, 571– Evans, E.; Ritchie, K. Dynamic Strength of 577. Molecular Adhesion Bonds. Biophys. J. 1997, 72, 1541–1555. [17] Takahashi, H.; Rico, F.; Chipot, C.; Sheuring, S. α-Helix Unwinding as Force Buffer in Grubm¨ uller, H. Force Probe Molecular DySpectrins. ACS Nano 2018, 12, 2719–2727. namics Simulations. Methods Mol. Biol. 2005, 305, 493–515. [18] Pronk, S.; Pall, S.; Schulz, R.; Larsson, P.; Bjelkmar, P.; Apostolov, R.; Shirts, M.; Pierse, C. A.; Dudko, O. K. DistinguishSmith, J.; Kasson, P.; van, d. S. D.; ing Signatures of Multipathway ConformaHess, B.; Lindahl, E. GROMACS 4.5: a Hightional Transitions. Phys. Rev. Lett. 2017, 118, Throughput and Highly Parallel Open Source 088101. Molecular Simulation Toolkit. Bioinformatics 2013, 29, 845–854. Yu, X.; Leitner, D. M. Vibrational Energy

Transfer and Heat Conduction in a Protein. [19] Best, R. B.; Hummer, G. Optimized MolecJ. Phys. Chem. B 2003, 107, 1698–1707. ular Dynamics Force Fields Applied to the Helix-Coil Transition of Polypeptides. J. Phys. [9] Barends, T. R. et al. Direct Observation of UlChem. B 2009, 113, 9004–9015. trafast Collective Motions in CO Myoglobin upon Ligand Dissociation. Science 2015, 350, [20] Lindorff-Larsen, K.; Piana, S.; Palmo, K.; 445–450. Maragakis, P.; Klepeis, J. L.; Dror, R. O.; Shaw, D. E. Improved Side-Chain Torsion Po[10] Brinkmann, L. U. L.; Hub, J. S. Ultrafast tentials for the Amber ff99SB Protein Force Anisotropic Protein Quake Propagation afField. Proteins 2010, 78, 1950–1958. ter CO Photodissociation in Myoglobin. Proc. Natl. Acad. Sci. U.S.A. 2016, 113, 10565– [21] Jorgensen, W. L.; Chandrasekhar, J.; 10570. Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of Simple Potential Functions [11] Young, H. T.; Edwards, S. A.; Gräter, F. How for Simulating Liquid Water. J. Chem. Phys. Fast Does a Signal Propagate through Pro1983, 79, 926–935. teins? PLoS ONE 2013, 8, e64746. [12] Aponte-Santamar´ıa, C.; Brunken, J.; [22] Joung, I. S.; Cheatham III, T. E. Determination of Alkali and Halide Monovalent Gräter, F. Stress Propagation through Ion Parameters for Use in Explicitly Solvated Biological Lipid Bilayers in Silico. J. Am. Biomolecular Simulations. J. Phys. Chem. B Chem. Soc. 2017, 139, 13588–13591. 2008, 112, 9020–9041. [13] Bilaniuk, N.; Wong, G. S. K. Speed of Sound in Pure Water as a Function of Temperature. [23] Hess, B. P-LINCS: A Parallel Linear Constraint Solver for Molecular Simulation. J. J. Acoust. Soc. Am. 1993, 93, 1609–1612. Chem. Theory Comput. 2008, 4, 116–122. [14] Heynmann, B.; Grubm¨ uller, H. AN02/DNPHapten Unbinding Forces Studied by Molecu- [24] Berendsen, H. J. C.; van Gunsteren, W. F. lar Dynamics Atomic Force Microscopy SimuMolecular Liquids. NATO ASI Series 1984, lations. Chem. Phys. Lett. 1999, 303, 1–9. 135, 475–500. 9



[25] Bussi, G.; Donadio, D.; Parrinello, M. Canonical Sampling Through Velocity Rescaling. J. Chem. Phys. 2007, 126, 014101. [26] Parrinello, M.; Rahman, A. Polymorphic Transitions in Single Crystals: A New Molecular Dynamics Method. J. Appl. Phys. 1981, 52, 7182–7190. [27] Darden, T.; York, D.; Pedersen, L. Particle Mesh Ewald: An Nlog(N) Method for Ewald Sums in Large Systems. J. Chem. Phys. 1993, 98, 10089–10092. [28] Grum, V. L.; Li, D.; MacDonald, R. I.; Mondragón, A. Structures of Two Repeats of Spectrin Suggest Models of Flexibility. Cell 1999, 98, 523–535. [29] Daday, C.; Kolˇsek, K.; Gräter, F. The Mechano-Sensing Role of the Unique SH3 Insertion in Plakin Domains Revealed by Molecular Dynamics Simulations. Sci. Rep. 2017, 7, 11669. [30] Kabsch, W.; Sander, C. Dictionary of Protein Secondary Structure: Pattern Recognition of Hydrogen-Bonded and Geometrical Features. Biopolymers 1983, 22, 2577–2637. [31] Zink, M.; Grubm¨ uller, H. Mechanical Properties of the Icosahedral Shell of Southern Bean Mosaic Virus: A Molecular Dynamics Study. Biophys. J. 2009, 96, 1350–1363. [32] Kappel, C.; Dölker, N.; Kumar, R.; Zink, M.; Zachariae, U.; Grubm¨ uller, H. Universal Relaxation Governs the Nonequilibrium Elasticity of Biomolecules. Phys. Rev. Lett. 2012, 109, 118304.

10


Page 10 of 12

Page 11 of 12 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


TOC figure

11



Frauke Gräter is head of the research group "Molecular Biomechanics" at the Heidelberg Institute for Theoretical Studies (HITS) and Professor at Heidelberg University. After her PhD studies at the Max Planck Institute for Biophysical Chemistry, Goettingen (with Prof. Grubmüller), and the Shanghai Institute of Materia Medica, Shanghai, she moved for a postdoc with Prof. Berne to Columbia University, New York. Prior to joining HITS in 2009, she has been junior group leader at the MPG/CAS Partner Institute for Computational Biology, Shanghai. She investigates how proteins have been designed to specifically respond to mechanical forces in the cellular environment or as a biomaterial, e.g. in the process of blood coagulation, in tendons, or in spider silk fibers. To this end, her group uses and further develops various simulation techniques from the molecular to the mesoscopic scale. Frauke Gräter has been awarded the Ada Lovelace award for HPC in 2017.


Page 12 of 12

How Fast Is Too Fast in Force-Probe Molecular Dynamics Simulations

Recommend Documents