An Ising-Like model for protein mechanical unfolding

The mechanical unfolding of proteins is investigated by extending the Wako-Saito-Munoz-Eaton model, a simplified protein model with binary degrees of freedom, which has proved successful in describing the kinetics of protein folding. Such a model is generalized by including the effect of an external force, and its thermodynamics turns out to be exactly solvable. We consider two molecules, the 27th immunoglobulin domain of titin and protein PIN1. In the case of titin we determine equilibrium force-extension curves and study nonequilibrium phenomena in the frameworks of dynamic loading and force clamp protocols, verifying theoretical laws and finding the position of the kinetic barrier which hinders the unfolding of the molecule. The PIN1 molecule is used to check the possibility of computing the free energy landscape as a function of the molecule length by means of an extended form of the Jarzynski equality.

The mechanical unfolding of proteins is investigated by extending the Wako-Saitô-Muñoz-Eaton model, a simplified protein model with binary degrees of freedom, which has proved successful in describing the kinetics of protein folding. Such a model is generalized by including the effect of an external force, and its thermodynamics turns out to be exactly solvable. We consider two molecules, the 27th immunoglobulin domain of titin and protein PIN1. In the case of titin we determine equilibrium force-extension curves and study nonequilibrium phenomena in the frameworks of dynamic loading and force clamp protocols, verifying theoretical laws and finding the position of the kinetic barrier which hinders the unfolding of the molecule. The PIN1 molecule is used to check the possibility of computing the free energy landscape as a function of the molecule length by means of an extended form of the Jarzynski equality. Manipulation experiments on single biomolecules have greatly increased our knowledge of the structural properties of such molecules. In a typical experiment a controlled force is applied to one of the free ends of the molecule, and the induced elongation is measured. Such experimental techniques have been used to probe the structure of proteins [1,2] and nucleic acids [3,4]. According to the common interpretation the unfolding of a molecule being pulled from one of its ends is hindered by kinetic barriers associated with the strongest linkages which serve to stabilize the molecular structure. The breaking of a molecular bond can thus be viewed as the overcoming of a kinetic barrier. It has been argued [5] that the study of the kinetics of bond breaking under different loading rates can provide much information about the internal structure of molecules, and in particular allows one to measure the strength of the molecular bonds, and to associate to them a position along the molecular structure. In the case of simple molecules, such as RNA hairpins, it is easy to obtain information on the molecular structure by pulling experiments [3]. However, when one deals with large molecular structures, such as multi-domain proteins, the inference of structural characteristics from the unfolding kinetics can prove a difficult task. Therefore, the study of simple models for the unfolding of proteins, whose microscopic native structures are known a priori, is highly desirable: investigating the kinetics of such models can shed new light on the relations between the experimentally observed unfolding features and the molecular structures.
The experiments discussed above are usually performed in non-equilibrium conditions: because of technical limitations the pulling process is faster than the typical molecular relaxation time. The problem of irreversibility of unfolding processes can be avoided by using the remarkable equality introduced by Jarzynski [6], which allows one to measure the free energy difference between the folded and the unfolded state of a biomolecule [7]. By using an extended form of the Jarzynski equality (JE) the free energy landscape of simple models of biopolymers has been probed as a function of the molecular elongation [8,9,10]. Although this approach appears very promising, it still has to be tested on systems of increasing complexity.
Here we approach the mechanical unfolding problem by means of a suitable generalization of the Wako-Saitô-Muñoz-Eaton (WSME) protein folding model [11,12,13]. This is a simplified statistical mechanical model where a binary variable is associated to each peptide bond. The equilibrium thermodynamics has been solved exactly [14,15] and the model has been quite successful in describing the kinetics of protein folding [16,17,18,19] and has also found applications in different fields (see [20] and references therein).
In the present paper we first extend the WSME model by considering the effect of an external force, and show that the equilibrium properties of this new model can be computed exactly, similarly to the original WSME model. In order to mimic the mechanical unfolding of proteins, we use computer simulations and study the unfolding kinetics of our model, both in the cases of constant force and dynamic loading. Finally we probe the free energy landscape of a model protein, exploiting the JE.
The WSME model describes a protein of N + 1 aminoacids as a chain of N peptide bonds (connecting consecutive aminoacids) that can live in two states (native and unfolded) and can interact only if they are in contact in the native structure and if all bonds in the chain between them are native. To each bond is associated a binary variable m k , with values 0, 1 for unfolded and native state respectively. The effective Hamiltonian of the model reads (1) where k B is the Boltzmann constant, and T is the absolute temperature. The first term assigns an energy ǫ i j < 0 to the contact (defined as in [13,14]) between bonds i and j if this takes place in the native structure (∆ i j = 1 in this case and ∆ i j = 0 otherwise). The second term represents the entropic cost q i > 0 of ordering bond i. In Ref. [14] it is shown how to compute exactly the partition function Z = {m k } exp[−βH 0 ({m k })] and the corresponding thermal averages. Here and in the following, the quantity β indicates the inverse temperature β = 1/(k B T ).
In order to couple the protein to an external force we assume that a configuration {m k } of the model defines a sequence of native stretches, separated by unfolded peptide bonds (see inset of Fig. 1). Peptide bonds i and j delimit a native stretch if and only if m i = m j = 0 and m k = 1 for i < k < j. A native stretch is regarded as a rigid portion of the molecule, even under application of the external force, with an end-toend length l i j . If j = i + 1 the stretch reduces to a single aminoacid, which is also regarded as a rigid structure, with length l i,i+i . The values of the parameters l i j are taken from the native structure of the protein [21]. Boundary conditions are introduced through the dummy bonds m 0 = m N+1 = 0. Given the direction of the external force, we assume that a native stretch, or a single aminoacid delimited by two successive unfolded bonds, can only take two orientations, parallel or antiparallel to the force, so that they contribute ±l i j to the length of the molecule. Therefore, given a configuration {m k } of the peptide bonds, we introduce a variable σ i j for each rigid portion of the molecule, either a native stretch or an aminoacid delimited by two successive unfolded bonds, taking values +1 (rigid portion parallel to the force) or −1 (antiparallel). Thus the end-to-end length of the molecule, in the force direction, reads Let us define the Hamiltonian H as the sum of the interaction energy term, contained in the effective Hamiltonian H 0 (1), and of the term − f L, which takes into account the effect of the external force Note that the definition of molecular length (2) is such that the set of variables {σ i j } is dynamically defined by the bond configuration {m k }: for each configuration {m k }, we consider only those variables σ i j such that and, as a function of {m k }, has the same structure as the WSME model (see Eq. (1)). Therefore its equilibrium thermodynamics is exactly solvable, as discussed in Ref. [14]. In the case f = 0, H eff reduces to Eq. (1) (up to an additive constant), with q i = ln 2. Thus, the effective entropic terms can be viewed as resulting from microscopic orientational degrees of freedom of the native stretches.
In the following, the quantity ǫ will indicate the system energy scale. This quantity, together with the interaction energies ǫ i j , will be chosen as in [13,14,18], see also [21]. We also introduce the reduced temperature in terms of the energy scale: T = k B T/ǫ . The time scale will be indicated by t 0 .
Following Ref. [14], for any choice of the model parameters, one can calculate the equilibrium value of the order parameter m = 1/N k m k and of the molecule length, as defined by Eq. (2), for varying force and temperature. We first consider the titin immunoglobulin domain I27 (89 aminoacids, pdb code 1TIT), which has been widely studied both experimentally [1] and theoretically [22]. The energy scale for such a molecule is taken to be ǫ/k B = 43.1 K, while the melting temperature is T m ≃ 8.03 [21]. In figure 1, the root mean square length of the 1TIT molecule is plotted as a function of the external force f for different temperatures. The plateau appearing in Fig. 1, in the low temperature regime, corresponds to the overall alignment of the molecule in the native configuration (m ≃ 1) to the force direction. Having introduced the molecular length we have built a framework within which the mechanically induced protein unfolding can be simulated by applying an external force. In the following, two manipulation schemes will be considered: in the first one the molecule is manipulated in a "force-clamp", where a sudden force is applied to one of the molecule's free ends, in the second one a time-varying force is applied, so that the load on the molecule increases gradually. In order to study the molecular unfolding, we run Monte Carlo simulations with Metropolis kinetics using the Hamiltonian H({m k } , {σ i j }, f (t)). In the following the time scale t 0 will correspond to a single Monte Carlo step. In the force clamp manipulation experiments, the molecule unfolds after a given time τ u which fluctuates between one realization of the unfolding process and the other, due to the stochasticity of the unfolding process. Unfolding can be viewed as an activated process [23], whose kinetics is dominated by a characteristic energy barrier ∆E u placed at the value x u of the reaction coordinate: therefore, it is usu-ally assumed that the average unfolding time τ u follows an Arrhenius law τ u = ω −1 0 exp β (∆E u − f x u ) , where ω 0 is a characteristic attempt rate depending on the microscopic features of the system.
We simulate force clamp manipulations of the 1TIT molecule as follows. Starting from thermal equilibrium with f = 0, at t = 0 we apply a non-zero force f and thus measure the unfolding time τ u as the first passage time of the molecule length across the threshold value L u , defined as half the molecule maximal length, L u = 140 Å, see figure 1. In figure 2 the unfolding time τ u of the 1TIT molecule, averaged over 1000 independent unfolding trajectories, is plotted as a function of the applied force f , for three values of the temperature T = 4, 6, 8. From a fit of the data shown in Fig. (2) to the Arrhenius law, we find that the unfolding length is x u ≃ 3 Å, and is independent of the temperature, as expected (see caption of fig. (2) for the exact values and uncertainties). Note that this value is in good agreement with the experimentally measured value of the titin unfolding length x u = 2.5 Å, found in Ref. [1]. We now consider the case where a time-dependent force is applied to our model molecule and the unfolding time is sampled over 1000 independent trajectories. Here the force increases linearly with time, with a rate r, and thus the rupture force f * is given by f * = rτ u , where unfolding time is defined as in the case of the force clamp. In refs. [23] it has been argued that, if the energy barrier ∆E u is large (compared to the thermal energy k B T ) and rebinding is negligible, the typical unbinding force of a single molecular bond under dynamic loading is given by where τ 0 is the characteristic unfolding time at zero force, τ 0 = ω −1 0 exp β∆E u . In Fig. 3 the breaking force f * is plotted as a function of the pulling velocity, for the 1TIT molecule, for three values of the temperature. The value of the unfolding length obtained by fitting the data to Eq. (4) is x u ≃ 3 Å, and is independent of the temperature, as expected (see caption of fig. 3 for the exact values). Note that this value agrees with that found with the force clamp manipulation, and with the experimental value x u = 2.5 Å found in Ref. [1]. On the other hand, in recent works [24], it has been argued that the rupture force f * has a more complex where the exponent ν depends on the microscopic details of the energy landscape, and γ is the Euler-Mascheroni constant γ ≃ 0.577. This equation reduces to Eq. (4) in the case ν = 1 or in the limit ∆E u → ∞ [24]. In the inset of Fig. 3, we plot the fits of the rupture force data to the equation definingf * . Although the agreement of the data with the equation appears to be rather good, the statistical errors of the fit parameters are quite large, sincef * depends nonlinearly on the set of the unknown parameters. This issue will be addressed in a forthcoming paper. We now evaluate the effective free energy landscape as a function of the molecular length L of the extended model here discussed. Formally, the free energy function where x t is the system microscopic configuration at time t, L(x) is the macroscopic length corresponding to microscopic configuration x, W t is the thermodynamic work done on the system by the external potential, up to the time t, defined by W t = t 0 dt ′ ∂H/∂t ′ , and the average is over all the trajectories of fixed duration t. In order to recover the partition function Z(L) from eq. (5) we use the procedure introduced and discussed in ref. [8]. The investigation of the free energy landscape of the 1TIT molecule, by using Eq. (5), turned out to be a very difficult task. Indeed, the typical value of the work associated to the unfolding of this molecule is of the order of some hundreds of k B T , for the value of ǫ here used. Since one has to evaluate exp(−βW t ), in order to exploit Eq. (5), the JE cannot give a reliable estimate of the energy landscape of the 1TIT molecule. For a discussion of the range of applicability of the JE to microscopic systems see, e.g., [25]. Therefore we consider the PIN1, a smaller protein whose folding characteristics have already been studied with the WSME model [18] (pdb code 1I6C, 39 aminoacids, ǫ/k B = 44 K [21]). Its free energy landscape as a function of the molecule elongation L, as given by Eq. (5), is plotted in Fig. 4 for different velocities of the pulling protocol. It can be seen that the curves F(L) collapse onto the same curve, as the pulling velocity is decreased. This is a clear signature that the energy landscape is correctly reconstructed, and its best estimate is the collapse curve, as discussed in Ref. [10]. On the other hand, the model introduced here is simple enough to allow the exact computation of the function Z(L), and hence we can obtain the exact value of the function F(L): the agreement with the landscape evaluated by the pulling manipulations is found to be very good, see Fig. 4.
In conclusion, we have introduced and studied a model of proteins under external loading. The unfolding length of the titin model is found to be in good agreement with the experimental one. We believe that this result represents a remarkable validation of the model that we have introduced: it suggests that our model, although minimal, captures the basic mechanisms underlying the unfolding process of proteins.
For a small protein, the extended form of the Jarzynski equality gives an estimate of the free energy landscape which is in good agreement with the expected one. We believe that our model can be successfully used to study the interplay between the protein structures and the kinetics of unfolding and refolding under external loading. As an example, by using computer simulations and applying an external force, one can easily determine which are the main contacts stabilizing the molecular structures. Furthermore, our model is also suitable to study the thermal unfolding of the proteins, in order to compare the thermal and the mechanical unfolding paths.
Appendix to "An Ising-like model for protein mechanical unfolding" In this appendix we briefly review the definition of the parameters ǫ, ǫ i j , ∆ i j in the WSME model, and discuss the definition of the new parameters l i j which have been introduced in the main text.
The parameters ∆ i j and ǫ i j , appearing in eqs. (1)-(3) of the main text are chosen following Ref. [1], starting from the protein native structure, as given in the Protein Data Bank (pdb in the following, http://www.pdb.org/). An atomic contact is present (∆ i j = 1) if, in the native state of the protein, at least two atoms from residues i and j + 1 (with j + 1 > i + 2) are closer than 4 Å. In this case ǫ i j is taken to be equal to kǫ, where k is an integer such that 5(k − 1) < n at ≤ 5k, and n at is the number of atomic contacts. As an example, in table I, the values of the contact parameter k, defined as k = ǫ i j /ǫ are listed for the 1TIT molecule.
The lengths l i j , in a generic N + 1 aminoacid protein, are defined as follows. Let us represent the aminoacid i with its N i − C α,i − C i sequence. Taking the native state as the reference configuration, l i j is chosen as the distance between the midpoint of the C i and N i+1 atoms and the midpoint of the C j and N j+1 atoms, see figure 5. If j = i + 1, l i,i+i corresponds to the length of a single aminoacid. If i = 0 the first point is substituted with the position of the N 1 atom, while when j = N + 1 we take the position of the C N+1 atom.
In order to fix the energy scale ǫ, we define the dimensionless unfolding temperature of our model T m , at zero force, as that temperature where the molecule order parameter m = 1/N k m k is equal to 1/2. For a given molecule, using the experimental value of the melting temperature T m , defined as the temperature where half of the sample is unfolded, we define the energy scale as ǫ ≡ T m k B / T m . This approach was used in Refs. [2,3].
Here and in the main text we indicate the proteins either with their common name or with their pdb code.
We find that the midpoint temperature of the 1TIT model molecule, at zero force, takes the value T ≃ 8.03. The experimentally measured unfolding temperature for such a protein is T m ≃ 346 K [4]. Consequently, we choose the energy scale to be ǫ/k B = 43.1 K.
The midpoint temperature of the PIN1 molecule (pdb code 1I6C), at zero force, is found to be T ≃ 7.55. The experimentally measured unfolding temperature for such a protein is T m ≃ 332 K [5]. Consequently, we choose the energy scale to be ǫ/k B = 44 K.