Dynamics of Disulfide-Bond Disruption and Formation in the Thermal

Sep 25, 2017 - ... by means of a series of coarse-grained canonical molecular-dynamics simulations, run at various temperatures, and replica-exchange ...
0 downloads 10 Views 1MB Size
Subscriber access provided by Purdue University Libraries

Article

Dynamics of Disulfide-Bond Disruption and Formation in the Thermal Unfolding of Ribonuclease A Pawel Krupa, Adam K. Sieradzan, Magdalena Anna Mozolewska, Huiyu Li, Adam Liwo, and Harold Abraham Scheraga J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.7b00724 • Publication Date (Web): 25 Sep 2017 Downloaded from http://pubs.acs.org on September 26, 2017

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

Journal of Chemical Theory and Computation 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 37

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

Journal of Chemical Theory and Computation

REVISED VERSION Dynamics of Disulfide-Bond Disruption and Formation in the Thermal Unfolding of Ribonuclease A

Paweł Krupa,1 Adam K. Sieradzan,2 Magdalena A. Mozolewska,1 Huiyu Li,1 Adam Liwo,2 and Harold A. Scheraga1*

1

Baker Laboratory of Chemistry and Chemical Biology, Cornell University, Ithaca, N.Y., 14853-1301,

U.S.A. 2

Faculty of Chemistry, University of Gdańsk, Wita Stwosza 63, 80-308 Gdańsk, Poland

*Corresponding author; phone: (607) 255 4034, fax: (607) 255 4700, e-mail: [email protected] 1 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

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

Page 2 of 37

Abstract Ribonuclease A (RNase A) is the most studied member of ribonucleases – a group of enzymes responsible for catalyzing RNA degradation. RNase A contains four disulfide bonds, which were found to be necessary for the native structure of the protein to form. In this work the kinetics and thermodynamics of RNase unfolding was studied by means of a series of coarse-grained canonical molecular-dynamics simulations, run at various temperatures, and replica-exchange molecular dynamics simulations with the UNRES force field, which is capable of dynamic formation and breaking the disulfide bonds during the course of the simulations. It was found that the Cys40-Cys95 bond was the first to break, while the Cys26-Cys84 bond was the last to break, independent of temperature, in agreement with available experimental data. Except for the Cys40-Cys95 bond, all disulfide bonds were relatively stable during the simulations. The formation/disruption of disulfide bonds was found to be temperature dependent for three out of four disulfide bonds in RNase A, except for the most stable disulfide bond between Cys65 and Cys72. A stable intermediate without the Cys40Cys95 disulfide bond, with structure similar to that of the most common folding intermediate observed experimentally, was found in simulated unfolding. In agreement with experiment, non- native disulfide bonds were also observed. By analyzing residue-position fluctuations, it was found that native disulfide bonds are located in the highly flexible regions of the protein, which is probably why their presence is necessary for the stability of RNase A.

Keywords: order of disulfide bond breaking, protein stability, molecular dynamics, UNRES force field 2 ACS Paragon Plus Environment

Page 3 of 37

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

Journal of Chemical Theory and Computation

1 Introduction Ribonucleases are widely spread enzymes in organisms, and hydrolyze ribonucleic acids. Bovine pancreatic ribonuclease A (RNase A) is perhaps the most extensively studied enzyme in the past century.1 It was selected as an enzyme worth studying because extensive prior investigations had already laid the groundwork for more precise experiments to address questions of current interest, starting with Anfinsen’s2 fundamental studies which connected the amino-acid sequence of the protein with its structure. The applications of disulfide-bond chemistry to studies of protein folding and unfolding, structure, and stability are usually illustrated with RNase A.3–8 RNase A is a single-domain protein with two-subdomain portions which contain 124 amino-acid residues comprising 3 α-helices, 7 short βstrands, and loops connecting them together (Figure 1). Additionally, its native structure contains the following four disulfide bonds: Cys26–Cys84, Cys40–Cys95, Cys58–Cys110, and Cys65–Cys72, and two cis-peptide bonds before prolines: Pro93 and Pro114, which stabilize the native structure of RNase A,4,9 making the folding process of RNase A slower10 and more difficult to study.11 RNase A is a model disulfide-bonded protein and it does not fold to its native structure until all of its native disulfide bonds have been formed;11 consequently, our effort here is focused on the unfolding process. Thermal unfolding was studied experimentally by Burgess & Scheraga as far back as 1975,12,13 initially showing that the C-terminus is relatively stable even at high temperatures. A few years later, it was found that, with increasing temperature, intermediates occur more frequently14. The most populated intermediate structure was found to be the one missing the Cys40-Cys95 disulfide bond.15 However, despite extensive experimental and theoretical work on the structure of intermediates16,17 and pathways and kinetics of the thermal unfolding of RNase A,18–21 many unexplained phenomena still exist. For that reason, in this work, the physics-based UNRES force field, which provides the

3 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

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

Page 4 of 37

possibility to form and break disulfide bonds dynamically during simulations22 was applied to study the unfolding process of RNase A extensively. UNRES has been tested successfully for protein structure prediction in the ‘Critical Assessment of protein Structure Prediction’ (CASP) blind-prediction tests23,24, and also proved to be able to treat protein folding25 and unfolding26 pathways. The unfolding simulations of RNase A provide strong evidence that disulfide-bond formation/disruption is temperature-dependent except for the disulfide bond Cys65-Cys72. Experimental studies suggest that the entropy of formation for the Cys65-Cys72 disulfide bond is the lowest of the four disulfide bonds;27,28 therefore, this bond is expected to be least dependent on temperature during the unfolding process, which is in agreement with our simulations (see section 3.3). In general, performed simulations offer the view that formation/disruption of disulfide bonds influences the structure of RNase A during the unfolding progress and that distinct order of disulfide-bond breaking can be observed.

2 Materials and Methods 2.1 UNRES model of polypeptide chain. UNRES29–32 is a coarse-grained model of polypeptide chains in which two interaction sites per residue are defined, namely the united side-chain (SC) and the united peptide group (p). The positions of the Cα atoms and those of the side-chain geometrical centers are used to define the geometry (Figure 2). The UNRES energy function is expressed by eq. 1; it constitutes an analytical approximation to the potential of mean force of polypeptide chains in water:  ∑     =  ∑   +  ∑    + 

   +    ∑    +

   ∑     + ! "  ∑ !  ,

$ 

+ % ∑ % &  +

4 ACS Paragon Plus Environment

Page 5 of 37

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

Journal of Chemical Theory and Computation

 ∑  '( , ) * + %+! ∑ %+! ,  + - " - + - . - + "

"

.

.

/+ " /+ + /+ . /+ + -   ∑"12 ∑ - '0

"

"

.

.

1

∑5 555

* + ∑5 !3+_55 + (1)

where the U's are energy terms, each multiplied by an appropriate weight, w; θi denotes the ith virtualbond valence angle, γi denotes the ith virtual-bond dihedral angle, and the angles, αSCi and βSCi define the location of the ith side chain with respect to the backbone (Figure 2). The term   denotes the mean free energy of the hydrophobic (hydrophilic) interactions between the side chains, which implicitly contain the contributions from the interactions of the side chain with the solvent. The terms   represent the excluded-volume potential of the side chain–

peptide group interactions. The terms  and   are the peptide-group interaction potential for   the Lennard-Jones interaction energy between peptide-group centers and the average electrostatic energy between peptide-group dipoles,33 respectively. The terms  and ! correspond to the

torsional and double-torsional potentials, respectively. The terms % ,  and %+! are the virtualbond angle-bending terms, the side-chain rotamer and the virtual-bond-deformation terms, respectively. The terms - are correlation or multibody contributions from the coupling between backbone-local 1

and backbone-electrostatic interactions, while the terms /+ are correlation contributions involving m 1

consecutive peptide groups; they are, therefore, termed turn contributions. The - terms are the

recently introduced knowledge- and physics-based side-chain backbone correlation potentials.31,32

2.2 Dynamic disulfide treatment in UNRES The term !3+_55 is a bimodal potential that accounts for switching between the bonded and nonbonded state of pairs of cysteine residues in water, described in ref 22, which facilitates simulations of 5 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

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

Page 6 of 37

dynamic formation of disulfide bonds. With the use of the dynamic treatment of disulfide bonds with the UNRES force field, all selected cysteine residues are able to form and break disulfide bonds with other cysteines forming both native and non-native bonds. The dynamic disulfide potential is given by equation 2 (see equation 4 and its description in ref. 22 for more details):

9  : ; < ; 7 =>>?@ ;B ; >>C$DE ?@ ; >>;B C F  >>G B 

!3+_55 =





   HI HI HI HI DE ?@ ; ;B C$J>K = ?@ ;B ; C F B G 

8 7 LM LM 6  : ; ≥ ;

(2)

where   is a harmonic-like bonded potential that also accounts for the anisotropy of cysteine 

residues (for details see equation 3 and its description in ref. 22), ;



between the centers of cysteine side chains in the bonded state, O



is the equilibrium distance

is the depth of the minimum

corresponding to two bonded cysteines, Ht is the height of the energy barrier between the non-bonded and bonded states of two cysteine residues, Rij is the distance between the centers of the side chains of two cysteine residues in the transition state,  O

LM

corresponding to two non-bonded cysteine residues, ;

LM

is the depth of the energy minimum is the corresponding equilibrium distance,

and g is a smoothing function connecting the two types of (bonded and non-bonded) potentials, which is given by equation 3 (see equation 8 and its description in ref. 22 for more details): TU 

PQ, R, S = @

C @3 − 2 %UC %U TU

(3)

For illustration of !3+_55 , see Figure 3 in ref. 22.

6 ACS Paragon Plus Environment

Page 7 of 37

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

Journal of Chemical Theory and Computation

The terms corresponding to factors of order higher than 1 are additionally multiplied by the respective temperature factors, fn(T), defined by eq 4:30

+  =

YZ [\]^$\]^ ]

b b YZ {\]^ab def g$\]^ [b def ]} c

(4)

c

where T0 = 300 K. It was found that !3+_55 could lead to formation of unphysical trisulfide bonds (three bonded cysteines) and higher order bond formation. To prevent formation of those bonds, an additional

potential 555 was added, and is described briefly in the next section.

2.3 A potential to prevent the formation of "trisulfide" bonds To prevent formation of trisulfide bonds, the following potential was designed and implemented for inclusion in the UNRES force field:

555 = ∑ ∈ ∑ m $; ∈ ∑lm $;l∈ , U i i



j $%

i $i

k $-

+

j



k

U' i * $%' $i * $-

U  i



j $%

 $ i 

k $-



+ (5)

where C is a set of specified cysteine residues, rij, rik and rjk are the distances between the side-chain centers of cysteine residues ij, ik and jk, respectively, whereas a, b, c and d are parameters, which have been optimized to produce a strong repulsion if three cysteine residues are at a contact distance simultaneously, and a nearly-zero contribution for two cysteine residues at a contact distance and the third residue farther away. As a model protein with multiple disulfide bonds, a hydrophobic cysteinerich protein from soybean (pdb code: 1HYP; Supporting Figure S1) was used to optimize the parameters of this newly implemented potential. The final parameters used were: a=50000.0 [Å-2]; 7 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

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

Page 8 of 37

b=0.05 [Å-6]; c=0.01; d=125000 [kcal×mol-1]. With these parameters, the native structure of the training protein was obtained in unrestricted simulations, with all disulfide bonds formed and no "trisulfide" bond appearing during the simulations. The trisulfide term Usss performs as follows. When the distances between all cysteines are large, the (rαβ+rβχ)6 terms (with α, β, γ equal to i, j, or k) in the denominators of the terms of equation 5 take over and the whole penalty function is effectively equal to 0 (Supporting Figure S2). When one of the distances is small, as happens for disulfide-bonded cysteines, the (rαβ-rβχ)2 term in one of the components of the penalty function will prevail only if one of the bonded cysteines also has a small distance (characteristic of a disulfide bond) with the third cysteine in the triad, a situation that should be prevented, by introducing an energy barrier.

2.4 Coarse-grained simulations

A series of canonical MD simulations was run with the coarse-grained UNRES force field at five temperatures: 278, 288, 298, 308 and 323 K, respectively, with and without dynamic treatment of disulfide bonds (see subsection 2.2). Each simulation consisted of 100 trajectories of 1 million steps, corresponding to approximately 4.6 µs real time per trajectory because of time-scale elongation resulting from the absence of fast-moving degrees of freedom in the coarse-grained treatment.29,34 Simulations were started from each of two initially energy-minimized structures of RNase A (PDB code: 1KF5), with RMSD values from the native structure equal to 0.207 and 0.180 Å for 1 and 10 cycles of optimization, respectively, to minimize the influence of the starting structure (mainly positions of united side-chains) on the results. Minimization was performed to remove steric overlaps resulting from the conversion of the experimental all-atom structure to the coarse-grained model. The VTS (variable time step)34 algorithm was used to control the length of the time steps. Temperature was controlled by the Langevin thermostat with the viscosity of water scaled down by a factor of 0.01 to 8 ACS Paragon Plus Environment

Page 9 of 37

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

Journal of Chemical Theory and Computation

speed-up the calculations, as in our earlier work.30,35 Snapshots were saved every 100 steps, providing 10 000 states per trajectory, and 1 million states per simulation to achieve satisfactory statistics during analysis.

To analyze the thermal stability of RNase A, an additional series of MD simulations was performed with the multiplexed replica exchange (MREMD)36–38 variant. Each simulation consisted of 144 trajectories of 1 million steps, corresponding to approximately 4.6 µs per trajectory, in the temperature range from 250 to 500K. As in canonical MD simulations, the Langevin thermostat controlled the temperature and the VTS algorithm the length of the time step. Exchange of replicas was attempted every 10 000 steps. The last half millions steps were used for the Weighted Histogram Analysis Method (WHAM)30,39,40 to determine the probabilities of the resulting structures at given temperatures to assess the fractions of folded (native-like) structures.

2.5 Data analysis 2.5.1 Free energy landscape analysis To study the free energy landscape, the second half of each trajectory was binned with use of the radius of gyration (Rg) and the Root-mean-square deviation (RMSD) from the native structure. The bin size was 0.2 Å for both variables. Subsequently, the free energy (∆F) was calculated by following eq 6:

∆p+ = −qrs ∑v d t t

uwf u

(6)

where ∆Fn is the free energy at the nth bin, R is the universal gas constant, T is the absolute

temperature of the thermostat, In is the number of structures in the nth bin, and ∑y U2 xU is the total

9 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

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

Page 10 of 37

number of structures. The graphs representing the free energy landscape are shown in Figure S3 of the Supporting Information.

2.5.2 Study of fluctuations of Cα atoms The fluctuations of the Cα atoms were calculated using the root-mean-square-fluctuation (RMSF) measure, in which the fluctuation of the atom j is given by equation 7.

qz{p = | ∑+ 2'; , − ;}, * 

+



(7)

where n is the number of snapshots in all trajectories, ; , is the position of atom j in snapshot i, and ;}, is the position of atom j in the native structure. The differences of positions between the snapshots and the native structure were calculated after optimal superposition of the structures.

2.5.3 Study of disulfide bond stability } The distance between the side-chain centers of two cysteine residues is ,3535 = 4.5 Å on

average; two cysteine side chains were considered as disulfide-bonded when the distance between their } side-chain centers was less than ,3535 + 1 Å = 5.5 Å. The distances between the side-chain centers

of the cysteine residues, which can form disulfide bonds were calculated from the last 0.5 million steps of each trajectory. A disulfide bond was considered as stable if ,3535 was below 5.5 Å for at least

75% of the time. Additionally, an analysis of the kinetics of breaking disulfide bonds along the simulations was performed, excluding the first 1% of the simulation, which corresponded to the relaxation of the system. Then the resulting fractions of broken disulfide bonds were fitted to the decay

10 ACS Paragon Plus Environment

Page 11 of 37

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

Journal of Chemical Theory and Computation

ratio of first order given by equation 8, assuming that all disulfide bonds are present at the beginning of the unfolding process (t = 0).

[~] =



l€ $l

‚1 − „ 'l€ $l * … = [~]† 1 − „ l 

(8)

with ‡ = ‡% + ‡F [~† ] = l

(9)



€ $lF

(10)

where [A] is the fraction of broken disulfide bonds, and kb and kf are the rate constants of the breaking and formation of disulfide bonds, respectively, k is the cumulative rate constant (observed in

experiment or in simulations), and ‚~† …is the equilibrium fraction of structures with the disulfide bond in question broken; k and [Aeq] can be obtained by fitting the first-order kinetic equation 8 to the simulation data, from which kb and kf can, in turn, be calculated by using equations 9 and 10. Subsequently ∆F of disulfide bond formation was calculated from [~]† , with the use of equation 11:

ˆpF1 = −qrs

[‰]Š‹ [‰]Š‹

(11)

Finally, ∆Fform was fitted to a linear equation to obtain the energy (∆E) and entropy (∆S) of each disulfide-bond formation:

ˆpF1 = ˆŒ − ˆ{

(12)

Because the dependence of ∆F on temperature is almost linear, we assumed that ∆E and ∆S are approximately constant in the range of temperatures considered.

11 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

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

Page 12 of 37

The fitted parameters of the kinetic equations and the simulated and fitted A[t] curves are shown in Table S1 and Figure S4 of the Supporting Information, while the calculated enthalpy and entropy changes are shown in Table S2 of the Supporting Information, respectively.

3 Results and Discussion 3.1 Thermal stability of disulfide bonds in the unfolding process Analysis of the disulfide-bond stability (Figure 3) shows that all disulfide bonds, except for the Cys40Cys95 disulfide bond are relatively stable at all temperatures. As mentioned in subsection 2.5.3 of the Methods section, the criterion for breaking a disulfide bond is the increase of the distance between the side-chain centers of the corresponding cysteine residues above 5.5 Å. It should be noted that disulfide bonds are breaking and forming multiple times during a simulation. The Cys40-Cys95 disulfide bond is stable in less than 50% of trajectories even at T = 278 K, which is the lowest examined temperature – and drops to 17% at T = 323 K; therefore, it breaks as the first one in most of the trajectories. The differences in the order of breaking the Cys58-Cys110 and Cys65-Cys72 disulfide bonds with increase of the temperature are also reflected in the stability of the disulfide bonds (Figure 3); the stability of the Cys58-Cys110 disulfide bonds increases slightly with increase of the temperature, while the stability of the Cys65-Cys72 bond is not temperature-dependent and remains similar at all temperatures. Therefore, the Cys58-Cys110 disulfide bond breaks after the Cys65-Cys72 bond at higher temperatures. The stability of the Cys26-Cys84 disulfide bond is high for lower temperatures and drops at 323K. Although the stability of this bond is lower than that of the Cys58-Cys110 and Cys65-Cys72 bonds, this bond breaks slower (Table 1), because of the high stability of the structural fragments in which both cysteine residues are involved, namely the second α-helix for Cys26 and the fourth β-strand for Cys84. These secondary-structure elements remain stable and conserved in the intermediate structures 12 ACS Paragon Plus Environment

Page 13 of 37

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

Journal of Chemical Theory and Computation

(compare Figure 1 and Figure 4B) also because of the favorable interactions between the polar and charged residues in the surroundings of both cysteines (Asn24, Tyr25, Asn27, and Gln28 near Cys26, and Thr82, Asp83, Arg85, and Glu86 near Cys84). Moreover, our studies show that only the presence of a single Cys26-Cys84 disulfide bond greatly stabilizes the structure of RNase A during unfolding, observed by a significant increase in the fraction of folded structures (Supporting Figure S5). This effect was not observed for any other disulfide bond.

3.2 Analysis of stability of RNase A and structure in the unfolding process Analysis of the RMSD and radius of gyration (Rg) shows that the presence of disulfide bonds increases the stability of RNase A, by decreasing the RMSD up to 2 Å; however, the chain remains flexible, because no noticeable change of Rg was observed (Supporting Figure S3). To identify the regions of RNase A, which contribute most to the fluctuations of the chain, root-mean-square fluctuations (RMSF) were calculated from the second half of each simulation (Figure 5). As was expected, the most fluctuating part of RNase A is the N-terminal fragment, which consists of three unstructured residues and the first α-helix (Figure 1). High fluctuations of the C-terminus can be observed only at T = 323 K, because the C-terminus is formed by the last β-strand (β7), which interacts strongly with other βstrands (β1, β4, β5, and β6) (Figure 1). The most flexible part of RNase A at low temperatures is the fragment composed of residues 64-70, which forms an outer loop, not interacting with other residues. High flexibility of the 64-70 fragment is clearly visible at higher temperatures, in which the Cys65Cys72 disulfide bond is necessary to stabilize this part of the protein. The 30-40 and 85-95 fragments constitute other very flexible parts of the RNase A structure and both of them are stabilized by the Cys40-Cys95 disulfide bond (Figure 5). However, due to the weakly defined character of these fragments, the Cys40-Cys95 disulfide bond is very fragile in a reducing environment (environment capable to disrupt disulfide bonds) and often breaks forming intermediate structures. With increase of 13 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

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

Page 14 of 37

temperature, the Cys85-Cys95 fragment becomes the primary source of fluctuations in the structure of RNase A, which affects the whole structure. Such a structure (Figure 4B) resembles the intermediate state found in the most common folding pathway of RNase A observed by Li et al.15,28 Using peptide mapping, another intermediate structure, lacking the Cys65-Cys72 disulfide bond, was observed by Rothwarf et al.41, which is also present in molecular dynamics simulations (Figure 4D). However, in general, the Cys65-Cys72 disulfide bond was found to be very stable in the simulations performed in this work and was also determined experimentally to be the most populated disulfide bond in one- and two-disulfide-bond intermediates.42–44 The flexible region of residues 16-25, which is the loop between the first and second α-helix, is stabilized by the Cys26-Cys84 disulfide bond at all investigated temperatures, keeping fluctuations on the same level in the 278-308 K temperature range. Also, the last loop in RNase A (residues 111-116) is greatly stabilized by the presence of the Cys58-Cys110 disulfide bond up to T = 308 K. Most of the intermediate structures preserve the secondary structure elements (Figure 4A-D); however, lack of one of the disulfide bonds can cause the separation of the RNase A subdomains, which is observed mostly for the structures that lack one disulfide bond, namely Cys65-Cys72 (Figure 4D). Such behavior has also been observed experimentally in the form of intramolecular distance distributions (IDDs), which are close to the experimental structure, but are significantly broader for the structures with one and more broken disulfide bonds.45,46 Higher flexibility of the N-terminal part of RNase A caused by looser packing and lack of disulfide bonds in the 1-25 residue fragment was also found experimentally.13,47

14 ACS Paragon Plus Environment

Page 15 of 37

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

Journal of Chemical Theory and Computation

3.3 Order of breaking disulfide bonds during thermal unfolding During the simulated thermal unfolding of RNase A a specific order of disrupting disulfide bonds was observed. The presence of multiple distinct unfolding pathways,15,48 and hence disulfide-bond disruption/formation is strictly connected to (un)folding pathways;49 a distinctive order of disulfidebond breaking, was also observed experimentally and theoretically for RNase A and other proteins48,50,51. At all temperatures, the first disulfide bond to break is the Cys40-Cys95 bond (in at least 80% of the trajectories; see Table 1), forming the intermediate structures, also observed experimentally as the most populated intermediate structures,15,28,52 especially in RNase A mutants,53 which can form native-like structures.54 Our kinetic analysis (for more details see section S1 ‘Kinetics and thermodynamics of disulfide bonds’ in Supporting materials) shows that the Cys40-Cys95 bond is, except for the simulation carried out at T = 278K, at least an order of magnitude faster to break (Supporting Table S1 and Supporting Figure S4). Conversely, Cys26-Cys84 was the last disulfide bond to break at all temperatures. Interestingly, even though this is not the most stable disulfide bond (Figure 3), our results suggest that it requires a cooperation of other structural elements to be disordered before this bond is broken (Table 1). Moreover, the Cys26-Cys84 bond has the most stabilizing effect of all disulfide-bonds in one-disulfide-bond mutants of RNase A (Supporting Figure S5), which could explain why disruption of the Cys26-Cys84 bond is the slowest. The Cys65-Cys72 bond is the most stable disulfide bond at low temperature while, at higher temperatures, the Cys58-Cys110 bond becomes the most stable one (Figure 3); therefore, the order of breaking the Cys58-Cys110 and Cys65Cys72 disulfide bonds depends on temperature. At T = 278 and 288K, the Cys58-Cys110 bond breaks before the Cys65-Cys72 bond for most of the trajectories (Table 1). With increasing temperature, the ratio of breaking the Cys58-Cys110 and Cys65-Cys72 bonds was close to 1 for T = 298 and T = 308K, while for T = 323K the Cys65-Cys72 bond broke before the Cys58-Cys110 bond (Table 1). RNase A 15 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

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

Page 16 of 37

with the Cys65-Cys72 disulfide bond broken was detected experimentally as the second most populated intermediate structure in thermal unfolding of this protein in a reducing environment.52 RNase A with the Cys40-Cys95 or Cys65-Cys72 disulfide bond broken was observed experimentally at wide ranges of temperatures and salt concentrations, while RNase A with the Cys26-Cys84 or Cys58Cys110 disulfide bond broken was stable only at lower temperatures (≤ 278 K).52

3.4 Analysis of MREMD simulations and non-native disulfide bonds An analysis of the MREMD simulations confirms that the presence of disulfide bonds stabilizes the structure of RNase A (Figure 6). In the simulations with dynamics treatment of the disulfide bonds, the peak of the fraction of folded structures increased over three times (from 0.12 to 0.38). With disulfide bonds, RNase A is very stable at the 270-295 K temperature range with the fraction of folded structures over 0.3. However, the structure is reasonably stable (fraction of folded structures above 0.15) even up to T = 330 K provided that the disulfide bonds are present (Figure 6). The results obtained here are in agreement with experimentally measured transitions temperatures, equal to 335 K for native RNase A and ~ 315 K for mutants lacking one disulfide bond at pH 6.4.15 Owing to the dynamic disulfide bond treatment in the simulation, which enables the cysteine residues to freely form bonds between all cysteine residues, both native and non-native, for some trajectories, some transient non-native disulfide bonds were formed. The most frequent was the Cys72Cys110 disulfide bond, then the Cys84-Cys95 bond and the Cys58-Cys72 bond, followed by the Cys26-Cys95, Cys58-Cys65 and, rarely, the Cys26-Cys40 bond. Non-native disulfide bonds are often found experimentally in the intermediate states in both protein folding and unfolding pathways and their formation and breaking seems to be necessary to achieve the native structure of a protein.49 Moreover, all of the above non-native disulfide bonds, except for the Cys26-Cys95 bond, have been 16 ACS Paragon Plus Environment

Page 17 of 37

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

Journal of Chemical Theory and Computation

observed by Vinci et al. in the early intermediates during experimental folding studies.55 However, Xu at el. observed these nonnative disulfide bonds experimentally in one-disulfide-bond intermediates at very high concentrations: Cys58-Cys65 (9.4%), Cys84-Cys95 (8.6%), Cys58-Cys72 (6.6%), Cys26Cys40 (4.0%), Cys72-Cys110 (2.0%), Cys26-Cys95 (1.1%).43 Reshuffling of disulfide bonds was also observed experimentally as one of the stages during unfolding for various proteins.56

4 Summary and Conclusions In this study, we investigated the thermal unfolding of RNase A by means of coarse-grained molecular dynamics with the UNRES force field that has the capacity of switching between the disulfide-bonded and non-bonded state of pairs of cysteine residues. It was found that, while the three disulfide bonds, namely Cys26-Cys84, Cys58-Cys110, and Cys65-Cys72, are generally stable and conserved in intermediate structures, the fourth disulfide bond between Cys40 and Cys95 was very often disrupted. This behavior can be explained by the fact that the Cys40-Cys95 bond is the only disulfide bond in RNase A, which connects two weakly defined regions, a long loop between the second α-helix and the first β-strand with another long loop between the fourth and the fifth β-strand (Figure 1). Moreover, the Cys40-Cys95 disulfide bond, contrary to the Cys26-Cys84 and Cys58-Cys110 bonds, which are buried inside the protein, is relatively exposed on the surface and therefore is more likely to be disrupted.15,28 Also the high flexibility of that region can cause the disruption of the disulfide bond and bending back of the loops. In particular, because Cys40 is surrounded by positively charged residues (Lys37, Asp38, Arg39 and Lys41), while the surrounding of Cys95 is more hydrophobic (Tyr92, Pro93, Asn94, Ala96, and Tyr97), the respective fragments dissociate from each other without a disulfide bond holding them together. In addition to the various stabilities of the disulfide bonds, the distinct order of disulfide-bond disruption was observed. At all investigated temperatures, the least stable Cys40-Cys95 disulfide bond 17 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

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

Page 18 of 37

was also the first to break and Cys26-Cys84 was the last one, while the order of breaking of Cys65Cys72 and Cys58-Cys110 disulfide bonds changes with temperature. Thermodynamic analysis revealed that entropy is a predominating factor which determines the stabilities of disulfide bonds. The disulfide bridges present in RNase A exhibit different temperature-dependent behavior. The Cys26-Cys84 disulfide bridge is stable at the temperature range from T = 278 K to T = 308 K and destabilizes at T = 323 K. The Cys40-Cys95 bridge is barely stable at the lowest temperature (T = 278 K) and becomes completely unstable at higher temperatures. The stability of the Cys58-Cys110 disulfide bond increases with increase of the temperature, with a peak at T = 308 K. The stability of the Cys65-Cys72 disulfide bond does not depend on temperature as long as the temperature is not too high, and this bond remains relatively stable even at comparatively high temperatures. However, the main unfolding pathways, including formation of temporary non-native disulfide bonds, remain similar at all investigated temperatures, which is consistent with experimental observations.57 Performed studies show that all four disulfide bonds are located in the most flexible parts of RNase A and their presence, even in a reducing environment, greatly stabilizes the most dynamic regions consisting of residues 16-25, 64-70, 85-95 and 111-116. Fluctuations of these loops, especially at higher temperatures, cause the instability of the whole molecule; therefore, the presence of all disulfide bonds is necessary for RNase A to achieve full stability. The disulfide bond which has the most stabilizing effect on the structure is the Cys26-Cys84 bond (Supporting Figure S5), which links two most fluctuating regions composed of residues 16-25 and 85-95. The role of the disulfide bonds increases with increase of the temperature, and RNase A with disulfide bonds is relatively stable up to T = 308 K. At T = 323 K the protein becomes unstable due to high kinetic energy of the chain causing high fluctuations of the whole structure.

18 ACS Paragon Plus Environment

Page 19 of 37

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

Journal of Chemical Theory and Computation

5 Acknowledgments This work was supported by grants from the U.S. National Institutes of Health (GM-14312) and the U.S. National Science Foundation (MCB10-19767), by grant UMO-2012/06/A/ST4/00376 from the Polish National Science Center (NCN). This research was supported by an allocation of advanced computing

resources

provided

by

the

U.S.

National

Science

Foundation

(http://www.nics.tennessee.edu/), and by the U.S. National Science Foundation through TeraGrid resources provided by the Pittsburgh Supercomputing Center. Computational resources were also provided by (a) the Informatics Center of the Metropolitan Academic Network (TASK) in Gdańsk, (b) the Interdisciplinary Center of Mathematical and Computer Modeling (ICM) of the University of Warsaw, (c) the OU Supercomputing Center for Education & Research (OSCER) at the University of Oklahoma (OU) (d) the 952-processor Beowulf cluster at the Baker Laboratory of Chemistry, Cornell University, and (e) the 692-processor Beowulf cluster at the Faculty of Chemistry, University of Gdańsk.

Supporting Information Supporting Information includes short description of the kinetics and thermodynamics of disulfide bonds (S1), representation of the 1HYP training proteins (Figure S1), 3D plot of the newly implemented potential preventing forming tri-sulfide bonds (Figure S2), 3D free energy plots of RNase A unfolding simulations (Figure S3), plots of the disrupted disulfide bonds (Figure S4), plots of the fractions of folded trajectories with one disulfide bond (Figure S5), results of fitting of the fractions of disulfide bonds disrupted during simulations to the first-order kinetic equation (Table S1), and calculated values of enthalpy and entropy of disulfide bonds (Table S2). This information is available free of charge via the Internet at http://pubs.acs.org 19 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

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

Page 20 of 37

References (1)

Marshall, G. R.; Feng, J. A.; Kuster, D. J. Back to the Future: Ribonuclease A. Biopolymers 2008, 90 (3), 259–277.

(2)

Anfinsen, C. B. Principles That Govern the Folding of Protein Chains. Science. 1973, 181, 223– 230.

(3)

Udgaonkar, J. B.; Baldwin, R. L. NMR Evidence for an Early Framework Intermediate on the Folding Pathway of Ribonuclease A. Nature 1988, 335 (6192), 694–699.

(4)

Schultz, D. A.; Baldwin, R. L. Cis Proline Mutants of Ribonuclease A. I. Thermal Stability. Protein Sci. 1992, 1 (7), 910–916.

(5)

Rothwarf, D. M.; Scheraga, H. A. Regeneration of Bovine Pancreatic Ribonuclease A. 1. Steady-State Distribution. Biochemistry 1993, 32 (10), 2671–2679.

(6)

Rothwarf, D. M.; Scheraga, H. A. Regeneration of Bovine Pancreatic Ribonuclease A. 2. Kinetics of Regeneration. Biochemistry 1993, 32 (10), 2680–2689.

(7)

Rothwarf, D. M.; Scheraga, H. A. Regeneration of Bovine Pancreatic Ribonuclease A. 3. Dependence on the Nature of the Redox Reagent. Biochemistry 1993, 32 (10), 2690–2697.

(8)

Rothwarf, D. M.; Scheraga, H. A. Regeneration of Bovine Pancreatic Ribonuclease A. 4. Temperature Dependence of the Regeneration Rate. Biochemistry 1993, 32 (10), 2698–2703.

(9)

Cook, K. H.; Schmid, F. X.; Baldwin, R. L. Role of Proline Isomerization in Folding of Ribonuclease A at Low Temperatures. Proc. Natl. Acad. Sci. U. S. A. 1979, 76 (12), 6157–6161.

(10)

Creighton, T. E. Kinetic Study of Protein Unfolding and Refolding Using Urea Gradient Electrophoresis. J. Mol. Biol. 1980, 137 (1), 61–80.

(11)

Neira, J. L.; Rico, M. Folding Studies on Ribonuclease A, a Model Protein. Fold. Des. 1997, 2 20 ACS Paragon Plus Environment

Page 21 of 37

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

Journal of Chemical Theory and Computation

(1), R1–R11. (12)

Burgess, A. W.; Weinstein, L. I.; Gabel, D.; Scheraga, H. A. Immobilized Carboxypeptidase A as a Probe for Studying the Thermally Induced Unfolding of Bovine Pancreatic Ribonuclease. Biochemistry 1975, 14 (2), 197–200.

(13)

Burgess, A. W.; Scheraga, H. A. A Hypothesis for the Pathway of the Thermally-Induced Unfolding of Bovine Pancreatic Ribonuclease. J. Theor. Biol. 1975, 53 (2), 403–420.

(14)

Matheson, R. R.; Scheraga, H. A. Steps in the Pathway of the Thermal Unfolding of Ribonuclease A. A Nonspecific Photochemical Surface-Labeling Study. Biochemistry 1979, 18 (12), 2437–2445.

(15)

Li, Y. J.; Rothwarf, D. M.; Scheraga, H. A. Mechanism of Reductive Protein Unfolding. Nat. Struct. Biol. 1995, 2 (6), 489–494.

(16)

Talluri, S.; Rothwarf, D. M.; Scheraga, H. A. Structural Characterization of a Three-Disulfide Intermediate of Ribonuclease A Involved in Both the Folding and Unfolding Pathways. Biochemistry 1994, 33 (34), 10437–10449.

(17)

Buckler, D. R.; Haas, E.; Scheraga, H. A. Analysis of the Structure of Ribonuclease A in Native and Partially Denatured States by Time-Resolved Non-Radiative Dynamic Excitation Energy Transfer between Site-Specific Extrinsic Probes. Biochemistry 1995, 34 (49), 15965–15978.

(18)

Zhang, J.; Peng, X.; Jonas, A.; Jonas, J. NMR Study of the Cold, Heat, and Pressure Unfolding of Ribonuclease A. Biochemistry 1995, 34 (27), 8631–8641.

(19)

Iwaoka, M.; Wedemeyer, W. J.; Scheraga, H. A. Conformational Unfolding Studies of ThreeDisulfide Mutants of Bovine Pancreatic Ribonuclease A and the Coupling of Proline Isomerization to Disulfide Redox Reactions. Biochemistry 1999, 38 (9), 2805–2815.

21 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

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

(20)

Page 22 of 37

Torrent, J.; Connelly, J. P.; Coll, M. G.; Ribó, M.; Lange, R.; Vilanova, M. Pressure versus Heat-Induced Unfolding of Ribonuclease A: The Case of Hydrophobic Interactions within a Chain-Folding Initiation Site. Biochemistry 1999, 38 (48), 15952–15961.

(21)

Narayan, M.; Xu, G.; Ripoll, D. R.; Zhai, H.; Breuker, K.; Wanjalla, C.; Leung, H. J.; Navon, A.; Welker, E.; McLafferty, F. W.; Scheraga, H. A. Dissimilarity in the Reductive Unfolding Pathways of Two Ribonuclease Homologues. J. Mol. Biol. 2004, 338 (4), 795–809.

(22)

Chinchio, M.; Czaplewski, C.; Liwo, A.; Ołdziej, S.; Scheraga, H. A. Dynamic Formation and Breaking of Disulfide Bonds in Molecular Dynamics Simulations with the UNRES Force Field. J. Chem. Theory Comput. 2007, 3 (4), 1236–1248.

(23)

He, Y.; Mozolewska, M.; Krupa P.; Sieradzan, A.; Wirecki, T.; Liwo, A.; Kachlishvili K.; Rackovsky, S.; Jagiełła, D.; Ślusarz, R.; Czaplewski, C.; Ołdziej, S.; Scheraga, H. A. Lessons from Application of the UNRES Force Field to Predictions of Structures of CASP10 Targets. Proc. Natl. Acad. Sci. U. S. A. 2013, 110 (37), 14936–14941.

(24)

Krupa, P.; Mozolewska, M. A.; Wiśniewska, M.; Yin, Y.; He, Y.; Sieradzan, A. K.; Ganzynkowicz, R.; Lipska, A. G.; Karczyńska, A.; Ślusarz, M.; Ślusarz, R.; Giełdoń, A.; Czaplewski, C.; Jagieła, D.; Zaborowski, B.; Scheraga, H. A.; Liwo, A. Performance of ProteinStructure Predictions with the Physics-Based UNRES Force Field in CASP11. Bioinformatics 2016, 32 (21), 3270–3278.

(25)

Maisuradze, G. G. G.; Medina, J.; Kachlishvili, K.; Krupa, P.; Mozolewska, M. A.; MartinMalpartida, P.; Maisuradze, L.; Macias, M. J. M. J.; Scheraga, H. A. H. A. Preventing Fibril Formation of a Protein by Selective Mutation. Proc. Natl. Acad. Sci. U. S. A. 2015, 112 (44), 13549–13554.

(26)

Tuszynska, I.; Bujnicki, J. M. Predicting Atomic Details of the Unfolding Pathway for YibK, a 22 ACS Paragon Plus Environment

Page 23 of 37

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

Journal of Chemical Theory and Computation

Knotted Protein from the SPOUT Superfamily. J. Biomol. Struct. Dyn. 2010, 27 (4), 511–520. (27)

Montelione, G. T.; Scheraga, H. A. Formation of Local Structures in Protein Folding. Acc. Chem. Res. 1989, 22 (2), 70–76.

(28)

Scheraga, H. A. Ribonucleases as Models for Understanding Protein Folding. In Ribonucleases; Nicholson, A. W., Ed.; Nucleic Acids and Molecular Biology; Springer Berlin Heidelberg: Berlin, Heidelberg, 2011; pp 367–397.

(29)

Ołdziej, S.; Czaplewski, C.; Liwo, A.; Chinchio, M.; Nanias, M.; Vila, J. A.; Khalili, M.; Arnautova, Y. A.; Jagielska, A.; Makowski, M.; Schafroth, H. D.; Kaźmierkiewicz, R.; Ripoll, D. R.; Pillardy, J.; Saunders, J. A.; Kang, Y. K.; Gibson, K. D.; Scheraga, H. A. Physics-Based Protein-Structure Prediction Using a Hierarchical Protocol Based on the UNRES Force Field: Assessment in Two Blind Tests. Proc. Natl. Acad. Sci. U. S. A. 2005, 102 (21), 7547–7552.

(30)

Liwo, A.; Khalili, M.; Czaplewski, C.; Kalinowski, S.; Ołdziej, S.; Wachucik, K.; Scheraga, H. A. Modification and Optimization of the United-Residue (UNRES) Potential Energy Function for Canonical Simulations. I. Temperature Dependence of the Effective Energy Function and Tests of the Optimization Method with Single Training Proteins. J. Phys. Chem. B 2007, 111 (1), 260–285.

(31)

Krupa, P.; Sieradzan, A. K.; Rackovsky, S.; Baranowski, M.; Ołdziej, S.; Scheraga, H. A.; Liwo, A.; Czaplewski, C. Improvement of the Treatment of Loop Structures in the UNRES Force Field by Inclusion of Coupling between Backbone- and Side-Chain-Local Conformational States. J. Chem. Theory Comput. 2013, 9 (10), 4620–4632.

(32)

Sieradzan, A. K.; Krupa, P.; Scheraga, H. A.; Liwo, A.; Czaplewski, C. Physics-Based Potentials for the Coupling between Backbone- and Side-Chain-Local Conformational States in the United Residue (UNRES) Force Field for Protein Simulations. J. Chem. Theory Comput. 2015, 11 (2), 23 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

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

Page 24 of 37

817–831. (33)

Liwo, A.; Czaplewski, C.; Pillardy, J.; Scheraga, H. Cumulant-Based Expressions for the Multibody Terms for the Correlation between Local and Electrostatic Interactions in the UnitedResidue Force Field. J. Chem. Phys. 2001, 115, 2323–2347.

(34)

Khalili, M.; Liwo, A.; Rakowski, F.; Grochowski, P.; Scheraga, H. A. Molecular Dynamics with the United-Residue Model of Polypeptide Chains. I. Lagrange Equations of Motion and Tests of Numerical Stability in the Microcanonical Mode. J. Phys. Chem. B 2005, 109 (28), 13785– 13797.

(35)

Khalili, M.; Liwo, A.; Jagielska, A.; Scheraga, H. A. Molecular Dynamics with the UnitedResidue Model of Polypeptide Chains. II. Langevin and Berendsen-Bath Dynamics and Tests on Model Alpha-Helical Systems. J. Phys. Chem. B 2005, 109 (28), 13798–13810.

(36)

Hansmann, U. Parallel Tempering Algorithm for Conformational Studies of Biological Molecules. Chem. Phys. Lett. 1997, 281, 140–150.

(37)

Mitsutake, A.; Sugita, Y.; Okamoto, Y. Replica-Exchange Multicanonical and Multicanonical Replica-Exchange Monte Carlo Simulations of Peptides. I. Formulation and Benchmark Test. J. Chem. Phys. 2003, 118 (14), 6664.

(38)

Czaplewski, C.; Kalinowski, S.; Liwo, A.; Scheraga, H. A. Application of Multiplexed Replica Exchange Molecular Dynamics to the UNRES Force Field: Tests with Alpha and Alpha+beta Proteins. J. Chem. Theory Comput. 2009, 5 (3), 627–640.

(39)

Kumar, S.; Rosenberg, J. M.; Bouzida, D.; Swendsen, R. H.; Kollman, P. A. THE Weighted Histogram Analysis Method for Free-Energy Calculations on Biomolecules. I. The Method. J. Comput. Chem. 1992, 13 (8), 1011–1021.

(40)

Gallicchio, E.; Andrec, M.; Felts, A. K.; Levy, R. M. Temperature Weighted Histogram 24 ACS Paragon Plus Environment

Page 25 of 37

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

Journal of Chemical Theory and Computation

Analysis Method, Replica Exchange, and Transition Paths. J. Phys. Chem. B 2005, 109 (14), 6722–6731. (41)

Rothwarf, D. M.; Li, Y. J.; Scheraga, H. A. Regeneration of Bovine Pancreatic Ribonuclease A: Identification of Two Nativelike Three-Disulfide Intermediates Involved in Separate Pathways. Biochemistry 1998, 37 (11), 3760–3766.

(42)

Talluri, S.; Falcomer, C. M.; Scheraga, H. A. Energetic and Structural Basis for the Preferential Formation of the Native Disulfide Loop Involving Cys-65 and Cys-72 in Synthetic Peptide Fragments Derived from the Sequence of Ribonuclease A. J. Am. Chem. Soc. 1993, 115 (8), 3041–3047.

(43)

Xu, X.; Rothwarf, D. M.; Scheraga, H. A. Nonrandom Distribution of the One-Disulfide Intermediates in the Regeneration of Ribonuclease A. Biochemistry 1996, 35 (20), 6406–6417.

(44)

Volles, M. J.; Xu, X.; Scheraga, H. A. Distribution of Disulfide Bonds in the Two-Disulfide Intermediates in the Regeneration of Bovine Pancreatic Ribonuclease A: Further Insights into the Folding Process. Biochemistry 1999, 38 (22), 7284–7293.

(45)

Navon, A.; Ittah, V.; Landsman, P.; Scheraga, H. A.; Haas, E. Distributions of Intramolecular Distances in the Reduced and Denatured States of Bovine Pancreatic Ribonuclease A. Folding Initiation Structures in the C-Terminal Portions of the Reduced Protein. Biochemistry 2001, 40 (1), 105–118.

(46)

Navon, A.; Ittah, V.; Scheraga, H. A.; Haas, E. Formation of the Hydrophobic Core of Ribonuclease A through Sequential Coordinated Conformational Transitions. Biochemistry 2002, 41 (48), 14225–14231.

(47)

McWherter, C. A.; Haas, E.; Leed, A. R.; Scheraga, H. A. Conformational Unfolding in the NTerminal Region of Ribonuclease A Detected by Nonradiative Energy Transfer. Biochemistry 25 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

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

Page 26 of 37

1986, 25 (8), 1951–1963. (48)

Juneja, J.; Udgaonkar, J. B. Characterization of the Unfolding of Ribonuclease A by a Pulsed Hydrogen Exchange Study: Evidence for Competing Pathways for Unfolding. Biochemistry 2002, 41 (8), 2641–2654.

(49)

Qin, M.; Wang, W.; Thirumalai, D. Protein Folding Guides Disulfide Bond Formation. Proc. Natl. Acad. Sci. U. S. A. 2015, 112 (36), 11241–11246.

(50)

Snijder, J.; van de Waterbeemd, M.; Glover, M. S.; Shi, L.; Clemmer, D. E.; Heck, A. J. R. Conformational Landscape and Pathway of Disulfide Bond Reduction of Human Alpha Defensin. Protein Sci. 2015, 24 (8), 1264–1271.

(51)

Zhang, L. Different Dynamics and Pathway of Disulfide Bonds Reduction of Two Human Defensins, a Molecular Dynamics Simulation Study. Proteins: Struct., Funct., Bioinf. 2017, 85 (4), 665–681.

(52)

Wedemeyer, W. J.; Welker, E.; Narayan, M.; Scheraga, H. A. Disulfide Bonds and Protein Folding. Biochemistry 2000, 39 (15), 4207–4216.

(53)

Xu, G.; Narayan, M.; Kurinov, I.; Ripoll, D. R.; Welker, E.; Khalili, M.; Ealick, S. E.; Scheraga, H. A. A Localized Specific Interaction Alters the Unfolding Pathways of Structural Homologues. J. Am. Chem. Soc. 2006, 128 (4), 1204–1213.

(54)

Xu, X.; Scheraga, H. A. Kinetic Folding Pathway of a Three-Disulfide Mutant of Bovine Pancreatic Ribonuclease A Missing the [40-95] Disulfide Bond. Biochemistry 1998, 37 (20), 7561–7571.

(55)

Vinci, F.; Ruoppolo, M.; Pucci, P.; Freedman, R. B.; Marino, G. Early Intermediates in the PDIAssisted Folding of Ribonuclease A. Protein Sci. 2000, 9 (3), 525–535.

26 ACS Paragon Plus Environment

Page 27 of 37

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

(56)

Journal of Chemical Theory and Computation

Chang, J. Y. A Two-Stage Mechanism for the Reductive Unfolding of Disulfide-Containing Proteins. J. Biol. Chem. 1997, 272 (1), 69–75.

(57)

Yan, Y. B.; Jiang, B.; Zhang, R. Q.; Zhou, H. M. Two-Phase Unfolding Pathway of Ribonuclease A during Denaturation Induced by Dithiothreitol. Protein Sci. 2001, 10 (2), 321– 328.

27 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

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

Page 28 of 37

Figure Legends Figure 1. Cartoon representation of RNase A (PDB code: 1KF5) with the cysteine residues involved in disulfide bonds marked by blue sticks and balls. Symbols α and β indicate α-helices and β-sheets, respectively. Figure 2. The UNRES model of polypeptide chains. Dark grey circles represent united peptide groups (p), white circles represent the Cα atoms, which serve to define backbone geometry, and light grey ellipsoids represent side chains (SCs). The p’s are located half-way between two consecutive Cα atoms. The virtual-bond angles θ, the virtual-bond dihedral angles γ, and the angles α and β that define the location of a side chain with respect to the backbone are also shown. Figure 3. Bar plots of the percentage of the trajectories, for which the disulfide bonds were stable during unfolding for Cys26-Cys84 (yellow), Cys40-Cys95 (red), Cys58-Cys110 (slanted blue), and Cys65-Cys72 (checkered green). Figure 4. Cartoon representation of the intermediate structures of RNase A with one disulfide bond broken: Cys26-Cys84 (panel A), Cys40-Cys95 (panel B), Cys58-Cys110 (panel C), and Cys65-Cys72 (panel D) in rainbow colors from N-terminus (blue) to C-terminus (red); bonded cysteine residues are shown in dark grey ball-and-stick representation and labeled, while non-bonded are light grey. Figure 5. Plots of root-mean-square fluctuations (RMSF) of each Cα atom, representing the flexibility of each residue, for simulations without disulfide bonds (solid blue line) and with dynamic treatment of disulfide bonds (solid black line). For clarity, the parts of the secondary-structure elements of the native structure have been marked in the bottom part of the diagram: α-helices (red lines), β-sheets (green lines), cysteine residues involved in disulfide bridges (yellow spots with dashed black lines indicating disulfide bonds).

28 ACS Paragon Plus Environment

Page 29 of 37

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

Journal of Chemical Theory and Computation

Figure 6. Diagrams of the fractions of folded trajectories (RMSD