Quantum Approach to Fast Protein-Folding Time

TL;DR

This paper proposes a quantum walk framework for protein folding, showing that quantum strategies can significantly reduce folding times compared to classical methods, offering new insights into the Levinthal paradox.

Contribution

It introduces a quantum approach to model protein folding as a quantum walk, avoiding prior hypotheses and demonstrating faster folding times.

Findings

Quantum walk model yields shorter folding times.

Quantum approach outperforms classical random walks.

Framework provides new insights into protein folding dynamics.

Abstract

In the traditional random-conformational-search model, various hypotheses with a series of meta-stable intermediate states were often proposed to resolve the Levinthal paradox. Here we introduce a quantum strategy to formulate protein folding as a quantum walk on a definite graph, which provides us a general framework without making hypotheses. Evaluating it by the mean of first passage time, we find that the folding time via our quantum approach is much shorter than the one obtained via classical random walks. This idea is expected to evoke more insights for future studies.

Equations

$${\hat{H}}_0=-\sum _{a,b}{J}_{ab}\mid s_a⟩⟨s_b\mid ,$$

$${\mathcal{E}}_a^{[\nu ]}=\sum _{k<l}{E}_{q_k{q}_l}{\delta }_{|{\text{r}}_k^a-{\text{r}}_l^a|,1}(1-{\delta }_{|k-l|,1}),$$

$$V^{[\nu ]}=\sum _a{\mathcal{E}}_a^{[\nu ]}\mid s_a⟩⟨s_a\mid .$$

$$\frac{\mathrm{d}}{\mathrm{d}t}{p}_a(t)=\sum _b{K}_{ab}{p}_b(t),$$

$${\tilde{T}}_{ab}=T_{ab}-{\Gamma }_{ab}+{\Lambda }_{ab}.$$

$$\frac{\mathrm{d}}{\mathrm{d}t}\hat{\rho }=\frac{1}{i\hslash }[\hat{H},\hat{\rho }]+\mathcal{L}(\hat{\rho }),$$

$$\mathcal{L}(\hat{\rho})=\frac{\lambda}{2}\bigl{(}2L\hat{\rho}\,L^{\dagger}-\hat{\rho}\,L^{\dagger}L-L^{\dagger}L\hat{\rho}\bigr{)}$$

$$P_{a,b}(t)={\int }_0^t{F}_{a,b}(t^{\prime })P_{b,b}(t-t^{\prime })\mathrm{d}{t}^{\prime }.$$

$${\tau }_{\mathrm{fd}}=\frac{{\int }_0^{{\tau }_0}t{F}_{1,c}(t)\mathrm{d}t}{{\int }_0^{{\tau }_0}{F}_{1,c}(t)\mathrm{d}t}.$$

Full text

Understanding how proteins fold spontaneously into their native structures is both a fascinating and fundamental problem in interdisciplinary fields involving molecular biology, computer science, polymer physics as well as theoretical physics etc.. Since Harrington and Schellman discovered that protein-folding reactions are very fast and often reversible processesSchellman1956 , there has been progressively more investigations on protein folding in both aspects of theory and experiment. Levinthal Levinthal noted early in 1967 that a much larger folding time is inevitable if proteins are folded by sequentially sampling of all possible conformations. Thus the protein was assumed to fold through a series of meta-stable intermediate states and the random conformational search does not occur in the folding process. The questions about what are the energetics of folding and how does the denature cause unfolding motivates one to think that the protein folding proceeds energetically downhill and loses conformational entropy as it goes. Based on such a hypothesis, the free-energy landscape framework was one way to describe the protein folding Malla ; Wolynes ; Passway:Jackson , where the energy funnel landscape provided a first conjecture of how the folding begins and continues Wolynes2012 .

As we known, there have been substantial theoretical models with different simplifying assumptions, such as Ising-like model Ising1 ; Ising2 , foldon-dependent protein folding model Englander , diffusion-collision model Karplus ; Sali , and nucleation-condensation mechanism Thirumalai ; Fersht etc.. Theoretical models are useful for understanding the essentials of the complex self-assembly reaction of protein folding, but till now, they often rely on various hypotheses Wolynes2005 ; Shakhnovich2006 ; Wolynes2012 ; Dill2012 ; Thirumalai2013 . This often brings in certain difficulties in connecting analytical theory to experimental results because some hypotheses can not be easily put into a practical experimental measurement. As it introduces less hypotheses in comparison to those theoretical models, the atomistic simulations PianaS ; HenryR ; Snow were used to investigate the protein folding along with nowadays’ advances in computer science. Recently, a high-throughput protein design and characterization method was reported that allows one to systematically examine how sequence determines the folding and stability Rocklin2017 . However, quantitatively achieving the folding time and accurately understanding how the sequence determines the protein folding remain to be a key challenge.

Here we propose a quantum strategy to formulate protein folding as a quantum walk on a definite graph, which provides us a general scheme without artificial hypotheses. In terms of the first-passage probability, one can calculate the folding time as the mean of the first-passage time. The obtained folding time in terms of our quantum scheme is much shorter than the one obtained via classical random walks. This idea is expected to open a new avenue for investigating the protein folding theoretically, which may motivate a necessary step toward developing technology for protein engineering and designing protein-based nanodevices Munoz .

Theoretical consideration: We describe the protein structure by the frequently adopted lattice model GoModel ; HP:Dill ; LiTang1996 ; LiJi2004 , namely, a protein is regarded as a chain of non-own intersecting unit (usually referring an amino-acid residue) of a given length on the two-dimensional square lattice. For a protein with $n$ amino-acid residues, we can calculate the total number $N_{n}$ of distinct lattice conformations that distinguish various protein intermediate structures. For instance, we have $N_{4}=4$ and $N_{6}=22$. This provides us a set with $N_{n}$ objects, which we call structure set and denote it by $\mathscr{S}_{n}=\{s_{1},s_{2},\cdots,s_{N_{n}}\}$ hereafter.

In order to study the protein folding process, we propose a concept of one-step folding. On the basis of the lattice model, we can naturally define the one-step folding by one displacement of an amino acid in one of the lattice sites. This makes us to establish certain connections between distinct points in the set $\mathscr{S}_{n}$ and to have a connection graph $\mathscr{G}_{n}$. In other words, two structures are connected via one-step folding if their conformation differs in one site only. As a conceptual illustration, we plotted the structure set $\mathscr{S}_{4}$, the connection graphs $\mathscr{G}_{4}$ and $\mathscr{G}_{6}$ in Fig. 1 (the $\mathscr{S}_{6}$ in Fig. S1 in supplementary material). Such a graph $\mathscr{G}_{n}$ is described by the so-called adjacency matrix $\mathrm{Mat}(J_{ab})$ that characterizes a classical random walk Kampen1997 on the graph.

Folding as a quantum walk: Letting $\mid s_{a}\,\rangle$ denote the state of a protein structure in the shape of the $a$-th lattice conformation, we will have a quantum Hamiltonian in a $N_{n}$-dimensional Hilbert space,

[TABLE]

where $J_{ab}$ refers to the connection between different points in the structure set, i.e., $J_{ab}$ is nonzero only if the $a$-th protein structure $s_{a}$ can be transited into the $b$-th structure $s_{b}$ by a one-step folding. With these physics picture one can also investigate quantum walk Aharonov1993 ; Farhi1998 ; WangJB2014 on the aforementioned graph.

From the coarse grained point of view, the 20 amino acids are classified HP:Dill as hydrophobic and hydrophilic (it is also called polar) groups according to their contact interaction. As $H$ and $P$ represent the hydrophobic and polar amino acids conventionally, a sequence of $n$ amino acids can be labeled by $\vec{q}=(q_{1},q_{2},\cdots,q_{n})$ where $q_{k}$ with $k=1,2,\cdots,n$ refers to either $H$ or $P$. Thus there will be totally a set of $2^{n}$ possible sequences. Let us call the entire of the whole random sequences as the sequence set denoted by $\mathscr{Q}_{n}=\{{\small[\nu]}\mid\nu=1,2,\cdots,2^{n}\}$. For any definite sequence-$[\nu]$ specified by a $\vec{q}$, we can calculate the total contact energy LiTang1996 ; LiJi2004 for each structure in $\mathscr{S}_{n}$,

[TABLE]

where $a=1,2,\cdots,N_{n}$ labels different structures, $k$ and $l$ denote the successive labels of the amino-acid residues in the sequence (i.e., the order in the chain), while $\textbf{r}^{a}_{k(l)}$ stands for the coordinate position of the $k(l)$-th residue in the $a$-th structure and $q_{k(l)}$ refers to either $H$ or $P$. Here the notation of Kronecker delta is adopted, i.e., $\delta_{\alpha,\beta}=1$ if $\alpha=\beta$ and $\delta_{\alpha,\beta}=0$ if $\alpha\neq\beta$. It is widely believed that the native structure of a protein possesses the lowest free energy Anfinsen . This can be interpreted by the hydrophobic force that drives the protein to fold into a compact structure with as many hydrophobic residues inside as possible HP:Dill . Thus the $H$-$H$ contacts are more favourite in the lattice model HP:Dill ; HP:Onuchic ; HP:Shakhnovich ; HP:Wolynes , which can be characterized by choosing $E_{PP}=0$, $E_{HP}=-1$, and $E_{HH}=-2.3$ as adopted in Ref. LiTang1996, .

With the contact energy (2) for every structure, the potential term can be expressed as

[TABLE]

Thus the total Hamiltonian for a definite sequence-$[\nu]$ is given by $\hat{H}^{[\nu]}=\hat{H}_{0}+V^{[\nu]}$. Clearly, the kinetic term $\hat{H}_{0}$ is determined by the connection graph $\mathscr{G}_{n}$ merely while the potential term $V^{[\nu]}$ defined on the structure set $\mathscr{S}_{n}$ is related to the concrete sequence-$[\nu]$ under consideration. This means that we have a hierarchy of Hamiltonian $\{H^{[\nu]}\mid\nu=1,2,\cdots,2^{n}\}$ actually for a theoretical study of the protein folding problem.

Note that one may obtain the same contact energy $\mathcal{E}_{a}$ for several different sequences. In this case, the dynamical properties are the same although those sequences may differ. Such a dynamical degeneracy implies a partition within the sequence set $\mathscr{Q}_{n}$. There are totally 16 possible sequences in $\mathscr{Q}_{4}$ which is partitioned into three subsets, i.e., $\mathscr{Q}_{4}=\{Q_{1},Q_{2},Q_{3}\}$, thus there will be three situations in the discussion on the time evolution. For $n=6$, there are totally 64 possible sequences in $\mathscr{Q}_{6}$ which is partitioned into 45 subsets, i.e., $\mathscr{Q}_{6}=\{Q_{1},Q_{2},\cdots,Q_{45}\}$ (see Tables SI $\&$ SII in supplementary material).

Random walk with sticky vertices: As we known, the continuous time classical random walk Montroll1965 on a graph $\mathscr{G}_{n}$ is described by the time evolution of the probability distribution $p_{a}(t)$ that obeys the master equation

[TABLE]

where $K_{ab}=T_{ab}-\delta_{ab}$ with $T_{ab}$ being the probability-transition matrix. In the conventional classical random walk, the probability-transition matrix is determined by the adjacency matrix of an undirected graph, namely, $T_{ab}=J_{ab}/\mathrm{deg}(b)$ where $\mathrm{deg}(b)=\sum_{c}J_{cb}$ represents the degree of vertex-$b$ in the graph $\mathscr{G}_{n}$. However, we ought to reconsider the random walk if there are some “sticky” vertices in the graph. This corresponds to the case when we take account of the contact energy $\mathcal{E}_{a}$ in the protein conformations. Thus, the probability-transition matrix should be modified so that the strength hopping into differs from that hopping out of those sticky vertices. The modified transition matrix $\tilde{T}$ is given by

[TABLE]

Here $\Gamma_{ab}=\theta({\cal E}_{ab})\Omega_{ab}T_{ab}$, $\Omega_{ab}={\cal E}_{ab}^{2}/({\cal E}_{ab}^{2}+1)$, and $\Lambda_{ab}=\delta_{ab}\sum_{c}\Gamma_{cb}$ where a notation ${\cal E}_{ab}={\cal E}_{a}-{\cal E}_{b}$ is adopted for simplifying the expression. The newly added two terms in (5) together guarantee the probability conservation. Therefore, in the presence of sticky vertices, one needs to solve the master equation (4) with the modified $\tilde{K}=\tilde{T}-I$ in the discussion of classical random walks.

The quantum dynamics: To accomplish a quantum mechanical understanding, we take account of the energy dissipation caused by the medium in which the folding occurs. This is governed by the Lindblad equation Lindblad

[TABLE]

where

[TABLE]

reflects the effect of dissipation. Here $L$ and $L^{\dagger}$ is called the Lindblad operator which can be determined from the analyses of random walks in the presence of sticky vertices. The aforementioned off-diagonal part $\Gamma$ in (5) provides this operator i.e., $L^{\dagger}=\sum_{ab}\Gamma_{ab}\mid s_{a}\,\rangle\langle\,s_{b}\mid$. Actually, equation (7) presents a general expression, which becomes the traditional one in terms of Pauli matrices, $\mathcal{L}(\hat{\rho})=(2\sigma^{-}\hat{\rho}\,\sigma^{+}-\hat{\rho}\,\sigma^{+}\sigma^{-}-\sigma^{+}\sigma^{-}\hat{\rho})\gamma/2$ with $\gamma=\lambda\Omega^{2}$ for a two level system that can be regarded as the two-vertices graph with a sticky vertex.

We solve the density matrix $\hat{\rho}(t)$ from Eq. (6) with the initial condition $\hat{\rho}(0)=\mid\!s_{1}\rangle\langle s_{1}\!\mid$. Here $\mid s_{1}\rangle$ refers to the completely unfolded straight-line structure. To illustrate our theory intuitively, we start from the simplest model of $n=4$ where the protein-folding problem becomes a task to investigate the quantum walk on the graph $\mathscr{G}_{4}$. We solve the $\hat{\rho}(t)$ numerically for the three situations $Q_{1}$, $Q_{2}$ and $Q_{3}$ respectively. In the calculation, we set $\hbar$ and $J$ to be unity and take the time step as $\Delta t=0.02$. For the initial condition: $\rho_{11}=1$ and the other matrix elements vanish when $t=0$, we solve Eq. (6) by means of Runge-Kutta method and obtain the magnitude of $\rho_{ab}^{(1)}(t)$ at any later time, $t=j*\Delta t$ with $j=1,2,\cdots$. We plot the time dependence of the diagonal elements of the solved density matrix for the $Q_{3}$ case in Fig. 2(a) and the other cases in the supplementary material Fig. S3. Likewise, we solve the density matrix for another initial condition $\rho_{44}(0)=1$ again so that the first-passage probability can be determined later on. The population of the most compact structure $\mid s_{c}\,\rangle$ is evaluated by the diagonal element $\rho_{cc}(t)$. For instance, $c=4$ in $\mathscr{S}_{4}$, and $c=9,19,$ and $20$ in $\mathscr{S}_{6}$. We can see that the probability of the state referring to the most compact structure $\mid s_{4}\,\rangle$ increases much more rapidly in the quantum folding process Fig. 2(a) than in the classical process Fig. 2(d). Toward a genuine understanding, we further study the quantum walk on the graph $\mathscr{G}_{6}$ by solving the density matrix numerically one by one for the aforementioned forty-five situations.

The folding time: Now we are in the position to define the protein folding time which can be formulated with the help of the concept of the mean first-passage time FirstPassageTime1969 ; FirstPassageTime2001 ; FirstPassageTime2004 ; FirstPassageTime2007 ; FirstPassageTime2016 . The mean first-passage time from a starting state $\mid s_{a}\,\rangle$ to a target state $\mid s_{b}\,\rangle$ is given by $\int_{t=0}^{\tau_{0}}tF_{a,b}(t)\mathrm{d}t\,/\int_{t=0}^{\tau_{0}}F_{a,b}(t)\mathrm{d}t$ where $\tau_{0}$ represents the time period when the first-passage probability vanishes $F_{a,b}(\tau_{0})=0$ which really occurs for the aforementioned quantum walk. For example in Fig. 2(b), the solved first-passage probability $F_{1,4}(t)$ becomes negative after $t=1.7$. The first-passage probability $F_{a,b}(t)$ from a state $\mid s_{a}\,\rangle$ to another state $\mid s_{b}\,\rangle$ after $t$ time obeys the known convolution relation

[TABLE]

Here $P_{a,b}(t)$ denotes the probability of a state being the basis state $\mid b\,\rangle$ at time $t$ if starting from the state $\mid a\,\rangle$ at initial time $t=0$. Quantum mechanically, it is evaluated by the diagonal elements of the density matrix, i.e., $P_{a,b}(t)=\rho^{(a)}_{bb}(t)$ where $\rho^{(a)}_{bb}(t)=\langle\,b\mid\hat{\rho}^{(a)}(t)\mid b\,\rangle$ is solved from Eq. (6) with the initial condition $\hat{\rho}(0)=\mid a\,\rangle\langle\,a\mid$, while $P_{b,b}(t)=\rho^{(b)}_{bb}(t)$ is solved with another initial condition $\hat{\rho}(0)=\mid b\,\rangle\langle\,b\mid$. Here the superscripts are introduced to distinguish the solution from different initial conditions. In the classical case, $P_{a,b}$ and $P_{b,b}$ refer to the $p_{b}(t)$ solved from Eq. (4), respectively, with initial conditions $p_{c}(0)=\delta_{ac}$ and $p_{c}(0)=\delta_{bc}$.

As protein folding is the process that proteins achieve their native structure, the folding time is the case that the starting state is chosen as $\mid s_{1}\,\rangle$ and the target states are the most compact states. For example, they are $\mid s_{9}\,\rangle$, $\mid s_{19}\,\rangle$ or $\mid s_{20}\,\rangle$ for $n=6$. The formula for the calculation of the folding time is thus given by

[TABLE]

To calculate the folding time we need to solve the first-passage probability $F_{1,c}(t)$ as a function of $t$ from the convolution relation (8). As an illustration, we first consider the case of $n=4$. For the classical folding process, we plot in Fig. 2(e) the $p^{(1)}_{4}(t)$ and $p^{(4)}_{4}(t)$. With these two time-dependent functions, the first passage-probability $F_{1,4}(t)$ can be further solved from the convolution relation (8) by numerical iterations (see Fig. 2(e)). It is nonnegative and approaches to zero when $t$ goes to infinity. This can be understood without difficulty because the classical probability distribution changes monotonously and approaches to its steady solution at the infinity time. However, for a quantum walk the probability distribution oscillates in time. We can see that the solved density matrix shown in Fig. 2(a) and Fig. S3 oscillates in time. With this new characteristics in quantum walk, the value of the first-passage probability solved directly from (8) appears to be negative in certain time region (see Fig. 2(b)) that is unphysical.

The zero point of $F_{1,4}(\tau_{0})=0$ determines the upper limit of the integration in the formula (9). In the simplest model with 4 residues, the classical folding times $\tau^{\mathrm{c}}_{\mathrm{fd}}$ for the sequence subsets $Q_{1}$, $Q_{2}$ and $Q_{3}$ are $6.0602$, $6.0351$ and $6.0180$ respectively. Their corresponding quantum folding times are $1.3208$, $1.2182$ and $0.9670$ respectively. Clearly, the quantum folding is faster than the classical folding with about four to six times even for the simplest model. In the same way, we calculate the quantum folding time for the forty-five situations for the case with 6 residues (see Tables SIII, SIV $\&$ SV in supplementary material). One can see that the quantum folding is faster than the classical folding with almost ten to hundred times or more. The experimental observation FoldingTimeExp ever exhibited that the protein folding is much faster than the theoretical prediction based on a random conformation search process. To visualize more easily we plot the quantum folding times $\tau_{\mathrm{fd}}$ in Figs. 2 (g) to (i). As a comparison, we also plot the ratios of classical folding time $\tau^{\mathrm{c}}_{\mathrm{fd}}$ to the quantum folding time $\tau_{\mathrm{fd}}$ on the same panels. In those three histograms, the longest folding time takes place for the sequence subsets $Q_{13}$, $Q_{38}$ and $Q_{42}$ while the shortest folding time occurs for the sequence subset $Q_{45}$, $Q_{29}$ and $Q_{29}$. The largest ratios $\tau^{\mathrm{c}}_{\mathrm{fd}}/\tau_{\mathrm{fd}}$ occur for the subsets $Q_{33}$, $Q_{31}$ and $Q_{41}$ but the smallest ratios occur for $Q_{17}$, $Q_{38}$ and $Q_{10}$.

In the above, we proposed a self-contained general theory to investigate protein folding problem quantum mechanically. In terms of $H$-$P$ lattice model, one can always have a structure set $\mathscr{S}_{n}$ for an amino-acid chain of any given number $n$ of residues. With such a structure set, one can naturally define a connection graph $\mathscr{G}_{n}$ by means of our definition of one-step folding. Thus either a classical random walk or a quantum walk on the graph can be solved with standard procedures. The former implies a random conformational search while the latter involves in fact a parallel search due to the quantum mechanical coherence Christopher . The application of quantum walk has attracted more attentions XuPeng to study various contemporary topics in recent years, our present strategy may open a new avenue in the area of the application of quantum walks. We have known if proteins were folded by sequentially sampling of all possible conformations, the calculated folding time would be inevitably very large because there is a very large number of degrees of freedom in an unfolded polypeptide chains. We elucidated that the quantum evolution naturally helps us to understand a faster protein folding. In terms of the concept of first-passage probability, we can calculate the quantum protein folding time as the mean first-passage time. It is worthwhile to mention that the first-passage probability solved from the conventional convolution relation may take negative value in some time domain. This is very important for the application of the quantum approach to an investigation of protein folding time. According to our results for $n=4$ and $6$, the quantum folding time is much shorter than that obtained from classical random walk. The presented theory is expected to bring in new insight features of protein folding process.

The work is supported by National Key R & D Program of China, Grant No. 2017YFA0304304, and partially by the Fundamental Research Funds for the Central Universities.