Role of pulling direction in understanding the energy landscape of proteins

Recent experiments have shown that the energy landscapes governing the mechanical unfolding of proteins are different from the energy landscapes governing the thermal or chemical unfolding rief ; busta ; tskhovrebova ; busta1 ; carrion1 . Until recently it was believed that all residues of any given protein contribute equally to the thermodynamic stability of globular proteins, with the stability of the protein structure being provided by its hydrophobic core serrano ; eriksson . It was therefore expected that the mechanical unfolding of a protein should be insensitive to the direction of the applied force. Surprisingly, recent experiments have shown that a protein’s resistance to unfolding depends strongly on the pulling direction brockwell ; carrion ; li ; matouschek . Varying the pulling direction also gives Angstrom precise information about the structure of biomolecules dietz ; dietz1 .

Force extension curves obtained by single molecule experiments (unfolding of proteins, stretching of DNA) are often described by the worm like chain model (WLC model) or freely jointed chain model (FJC model) fixman ; doi . The mechanical properties of the FJC and WLC models are well understood, and the usefulness of these models stems from their simplicity, which means that the analytic expression for the extension can be written in a very simple form Marko ; rief98 . The only parameter appearing in the description of these models is the persistence length of the polymer chain. The force-extension curves obtained from these models Marko ; rief98 ; degennes show an excellent agreement with the main features observed in early experiments rief ; busta ; tskhovrebova ; busta1 . Since the FJC and WLC models are isotropic, in order to study the resistance of a protein to unfolding when the pulling geometry is changed, one would have to use different persistence lengths (along different directions) to fit the force extension curves brockwell ; carrion ; rief98 . So when a polymer has been pulled along the $x$-direction one may use a certain persistence length which reproduces the observed force-extension curve. If the pulling direction is changed another persistence length is used to fit the new force-extension curve. However, the underlying model is still isotropic and the procedure is somewhat ad hoc and does not provide a unified approach to the modeling of the anisotropic behavior creighton ; rosa seen in experiments. Equally importantly these simple models do not incorporate excluded volume effects in its description.

In a recent paper Kumar and Giri kumar used partially directed self-avoiding walks (PDSAWs) to model anisotropic biomolecules. When a force is applied to one end of the chain it undergoes a transition from a folded to an unfolded state. Numerical studies based on a chain size of $N=30$ showed that when the chain is being pulled along the preferred direction of the polymer, the nature of the transition is akin to unzipping, but when force is applied perpendicular to the preferred direction, the transition is akin to shearing.

In section $II$ we provide the analytical solution of the model and show that a change in the pulling direction gives rise to different phase diagrams even in the thermodynamic (infinite chain length) limit. In section  $III$ we report on an algorithmic breakthrough which has enabled us to increase the chain length to $N=150$ which is five times longer than in the previous study kumar . The longer series are used to obtain numerical estimates for the phase boundary in excellent agreement with the exact results. In this section we also discuss the relevance of two ensembles namely, the constant force ensemble (CFE) and the constant distance ensemble (CDE) which are appropriate when describing different experimental setups. With the enhanced information about the exact density of states of longer chains, we are able to study precisely finite-size effects which are crucial to all single molecule experiments. In particular we find marked differences between the cases where the force is applied along the $x$- or $y$-directions, respectively, and we explain these observed differences using simple heuristic arguments. The PDSAW is a very restricted model and in order to further examine the role of anisotropy we introduce an anisotropic self-attracting self-avoiding walk (ASASAW) model and using exhaustive exact enumeration data for chains up to length $N=48$ study the combined effect of anisotropy and pulling in section $IV$. The paper ends with a brief discussion in section $V$.

The PDSAW is a self-avoiding walk in which steps along some direction (say, the negative $x$-direction) are forbidden. We study this model on a two dimensional square lattice. An attractive energy $\epsilon_{u}$ is associated with each non-bonded nearest neighbor vander . At very low temperatures the attractive interactions dominate, and the PDSAW is in a collapsed phase, where the density of monomers in the bulk is close to one. The typical configurations of the chain mimics the $\beta$-sheet creighton (see Fig. 1). At high temperatures thermal fluctuations are strong enough to break some of the interaction bonds and with increasing entropic contributions to the free energy the chain can be in an extended phase. A force $\vec{F}$ may be applied along some fixed direction giving rise to a stretching energy $E_{s}$ arising from the applied force $\vec{F}$

where $\vec{\alpha}$ is the end-to-end vector. In this paper we consider only the cases where the force is applied either along the $x$-direction ($f_{x}$) or the $y$-direction ($f_{y}$).

The phase boundary between the collapsed and the extended phase can be obtained by calculating the macroscopic shape of the collapsed phase at low temperatures. Let $\omega_{x}=\exp(f_{x}/T)$, $\omega_{y}=\exp(f_{y}/T)$ and $u=\exp(\epsilon_{u}/T)$, where $T$ is the temperature. We set the Boltzmann constant $k_{B}$ equal to $1$. The energy of the configuration which mimics the $\beta$-sheet is (see Fig. 1)

where, as illustrated in Fig. 1, $a_{1}=|y_{1}-y_{0}|$ is the height (number of steps in the $y$-direction) of the first column, $a_{2}=|y_{b}-y_{b-1}|$ is the height of the last column and $b$ is the number of steps in the $x$-direction.

Following the approach developed in suppl , we express the sum appearing in Eq. (2) as a surface tension term. The macroscopic shape of the polymer is calculated and the value of $\omega_{x}$ at which the collapsed phase becomes unstable gives the critical value, $\omega_{x}^{c}$, in terms of $\omega_{y}$ and $u$:

Eq. (3) gives the complete three dimensional phase boundary. For $f_{y}=0$ ($\omega_{y}=1$) Eq. (3) reduces to the following simple expression for the critical force when pulling along the $x$-direction

This coincides with the result obtained by transfer matrix methods rosa . However, that method is not easily generalized to the case where force is applied along the $y$-direction. Eq. (3), on the other hand, is sufficiently general to yield an expression for the critical force when pulling along the $y$-direction. Setting $f_{x}=0$ (or $\omega_{x}=1$) yields the critical force when pulling along the $y$-axis:

For zero force (i.e. $\omega_{x}=\omega_{y}=1$) the expression reduces to the exact value $u^{c}=3.3829\ldots$ ($T=0.8205\ldots$ if we set $\epsilon_{u}=1$ ) at which thermal unfolding occurs. The exact phase boundaries obtained from Eqs. (4) and (5) are shown in Fig. 2 and compared to numerical results obtained for $N=150$.

Our numerical studies of the process of force induced protein unfolding are based on the exact enumeration of all possible conformations of the chain. This approach not only enables us to obtain precise estimates for the phase boundary, but also allows us to study the intermediate states crucial to events occurring during the mechanical stretching process altmann ; rief1 ; schwaiger . In this study we have combined the power of parallel processing and transfer matrix calculations to extend the limit of earlier studies not just by one or two monomers but five times over! Thus we have extended from the previous longest chain length kumar $N=30$ to much longer chains of up to $N=150$ steps. The partition function is

where $C_{N}$ is the total number of $N$-stepped PDSAWs having $N_{p}$ non-bonded nearest neighbor pairs with end-to-end vector $(\alpha_{x},\alpha_{y})$.

The value of the transition temperature (for fixed force and finite $N$) can be obtained from the fluctuations in the number of non-bonded nearest neighbor interactions. The fluctuations are defined as $\chi=\langle N_{p}^{2}\rangle-\langle N_{p}\rangle^{2}$, with the $k$’th moment given by

When the fluctuations $\chi$ are plotted as a function of temperature the resulting graph (Fig. 3a) has a peak at the transition temperature. The inset shows that the height of the peak in the fluctuation curve grows as a power-law with $N$, this being the hall mark of a phase transition. By setting $\epsilon_{u}=1$ we obtain (for $N=150$) at zero force ($\omega_{x}=\omega_{y}=1$) the transition temperature $T=0.8446$. This is already quite close to the value $T=0.8205$ obtained in the thermodynamic limit ($N\to\infty$). The bottom panel in Fig. 3 shows the variation in the transition temperature plotted against $1/N$ and we see a pronounced curvature for large $N$ making an extrapolation to the limit $N\to\infty$ very difficult.

In Fig. 2 we show the force-temperature diagrams for $N=150$ when pulling along the $x$- and $y$-axis (circles and squares, respectively). These numerical estimates for the phase-boundary are in excellent agreement with the analytical results (shown as solid lines) obtained from Eqs. (4) and (5). Thus exact enumeration technique works quite well even for small $N$ as far as the phase diagram is concerned. This has also been seen earlier mishra in case of the surface adsorption.

The model described above has an exact solution and gives qualitative features similar to experiments. However, the physical constraints imposed by experimental set-ups has not been taken fully into account. For example, in experiments using atomic force microscopes, receptor and ligand molecules are attached to a substrate and a transducer, respectively, which introduce anisotropy in the systems. Also steric constraints due to the confinement imposed by the experimental set-up can lead to a loss of entropy and thus result in modifications to the behavior of the chain. Furthermore the constraint of partially directing the walk so it cannot take steps in the negative $x$-direction is not really appropriate to most experimental conditions. Such a constraint is typically only appropriate when the system is under flow or constant external field lemak ; lemak1 . A more realistic model of polymers is provided by self-attracting-self avoiding walks (SASAWs), which is well suited to the modeling of a linear chain in a poor solvent degennes ; vander . However, this model is isotropic in nature. In this paper we introduce anisotropy to the model. We achieve this by considering nearest neighbor interactions with different strengths along the $x$- and $y$-directions as shown in Fig. 7a. This is in accordance with real proteins where the interactions along different directions can differ by orders of magnitude creighton . Hence by changing the strength of the interactions one can vary the degree of anisotropy in the system and study the effect this has on the force-temperature phase diagram. The partition function for anisotropic-self-attracting-self-avoiding walks (ASASAWs) may be written as

Here $C_{N}$ is the total number of ASASAWs of $N$ steps having $N_{px}$ and $N_{py}$ nearest neighbor pairs along the $x$- and $y$-directions, respectively, while $u_{x}=\exp[-\beta\epsilon_{x}]$ ($u_{y}=\exp[-\beta\epsilon_{y}]$) are the Boltzmann weights associated with the nearest neighbor interactions between non-bonded monomers along, respectively, the $x$- and $y$-directions.

The qualitative as well as quantitative features of the force-temperature phase diagram as shown in Fig. 8 for the isotropic case (i.e. $u_{x}=u_{y}=u$) remain the same as those reported in kumar1 . To make our model closely resemble the PDSAW discussed above, we set either $u_{x}=1$ or $u_{y}=1$ and apply a force along the $x$-direction. For $u_{y}=1$ the ground state is dominated by configurations similar to the one shown in Fig. 7(b). This closely resembles the PDSAW model (Fig. 1) subjected to a force along the $x$-axis and the nature of the transition is akin to an unzipping transition. On the other hand for $u_{x}=1$ the ground state is dominated by configurations like the one shown in Fig. 7(c). This closely resembles the configurations of Fig. 1 with a force along the $y$-axis and consequently the nature of the transition is akin to a shearing type transition. It should be noted that the number of ASASAWs, $C_{N}$, is different from the number of PDSAWs even if we put $u_{y}=1$. We observe a significant decrease in the transition temperature compared to SASAWs. This is caused by a decrease in entropy which allows the polymer to acquire conformations similar to a $\beta$-sheet. The effect of surface confinement is also evident from these plots. The transition temperature is found to be less for $u_{x}=1$ in comparison to $u_{y}=1$. Another major difference between the ASASAW and SASAW phase-diagrams is that there appears to be no re-entrance in the ASASAW model when $u_{x}=1$ or $u_{y}=1$. It is well known that the re-entrant behavior of the SASAW model is due to a non-zero ground state entropy, that is, the number of configurations in the ground state grow exponentially with $N$ (these configurations are just Hamiltonian walks inscribed in a rectangle with side lengths as close as possible to $\sqrt{N}$). For the ASASAWs we study here the ground-states are, as argued above, folded sheets as in Fig.7(b) and (c), respectively. The number of such configurations do not grow exponentially with $N$ and are thus too few to give a non-zero ground state entropy, which is why we do not see re-entrance. However, if both $u_{x}$ and $u_{y}$ are greater than 1 then one may see a re-entrant phase diagram.

The force extension curve shown in Fig. 9 not only shows the existence of intermediate states but also the emergence of new states depending on the direction of the interactions in accordance with the findings for PDSAWs. As the temperature is increased we observe in Fig. 10 (as expected) that the intermediate states are washed out due to the resulting increase in entropy. We also observe that the system attains a higher stability against the force when the interactions are along the $y$-direction and that intermediate states survive even at high temperatures as may be seen from Fig. 10.

Finally we analyze the ASASAW model in the constant distance ensemble. In Fig. 11 we have plotted the average force $\langle F\rangle$ kumar1 as a function of the extension when pulling in the $x$-direction with interactions either along the $x$- or $y$-direction. The results are qualitatively very similar to those obtained for the PDSAW model (see Fig. 6). In particular we note that the oscillation at low-temperature are much more pronounced in the case where the interactions are in the $y$-direction. As explained above (see Fig. 7) there is a close correspondence between the ground-state (low-temperature) configurations of the ASASAW with interactions in the $x$- or $y$-directions and the PDSAW with a force applied along the $x$- or $y$-direction, respectively. This close correspondence is reflected in the similarity of the results.

In summary we have demonstrated via the analytical solution for the phase boundary of PDSAWs as well as through high precision numerical calculations that changing the pulling direction can change the nature of the unfolding transitions and give rise to very different phase diagrams. The emergence of many new intermediate states suggests that there may be many pathways to the unfolding of a polymer. The findings for the PDSAW model have been confirmed by our study of ASASAWs. The inclusion of anisotropy in the traditional SASAW model of polymers gives unequivocal evidence that the features observed here for the PDSAW model are in fact generally true and not an artifact of the PDSAW model. In future additional work like extensive Monte Carlo simulations and molecular dynamics studies are needed to better understand the role played by anisotropy and how varying the pulling direction can help us gain a better and deeper understanding of biomolecules.

This research has been supported by the Department of Science and Technology and University Grants Commissions, India and also by the Australian Research Council. The calculations presented in this paper used the computational resources of the Australian Partnership for Advanced Computing (APAC) and the Victorian Partnership for Advanced Computing (VPAC).