Reaction Pathway Sampling and Free-Energy Analyses for Multimeric Protein Complex Disassembly by Employing Hybrid Configuration Bias Monte Carlo/Molecular Dynamics Simulation.

Physicochemical characterization of multimeric biomacromolecule assembly and disassembly processes is a milestone to understand the mechanisms for biological phenomena at the molecular level. Mass spectroscopy (MS) and structural bioinformatics (SB) approaches have become feasible to identify subcomplexes involved in assembly and disassembly, while they cannot provide atomic information sufficient for free-energy calculation to characterize transition mechanism between two different sets of subcomplexes. To combine observations derived from MS and SB approaches with conventional free-energy calculation protocols, we here designed a new reaction pathway sampling method by employing hybrid configuration bias Monte Carlo/molecular dynamics (hcbMC/MD) scheme and applied it to simulate the disassembly process of serum amyloid P component (SAP) pentamer. The results we obtained are consistent with those of the earlier MS and SB studies with respect to SAP subcomplex species and the initial stage of SAP disassembly processes. Furthermore, we observed a novel dissociation event, ring-opening reaction of SAP pentamer. Employing free-energy calculation combined with the hcbMC/MD reaction pathway trajectories, we moreover obtained experimentally testable observations on (1) reaction time of the ring-opening reaction and (2) importance of Asp42 and Lys117 for stable formation of SAP oligomer.

Introduction

The physicochemical entity of various biological phenomena is multimeric biomolecule complexes formed transiently or stably in the cell. The molecular components of such complexes have been extensively identified and individually assigned to the corresponding biological functions (e.g., signal transduction,[1−3] DNA repairs,[4,5] protein systhesis,[6,7] and RNA splicing[8,9]). Since these functions are expressed and regulated through biomolecule assembly and disassembly in response to the cellular environment,[10−12] physicochemical characterization of these processes is a milestone to understand the mechanisms for biological phenomena at the molecular level.

Multiple biomacromolecule assembly and disassembly processes can be physicochemically characterized by solving the two essential problems coupled with each other: (1) identifying a set of intermediate states characterized by subcomplexes appearing in these processes and (2) elucidating physicochemical mechanisms of transition between such intermediate states.

Mass spectrometry (MS) analyses have become a powerful tool to systematically analyze relative subunit interface strength and provide a possible set of subcomplexes undergoing assembly and disassembly processes at the molecular level.[13−15] Meanwhile, assembly and disassembly orders have been extensively studied in the field of structural bioinformatics (SB) during the last 15 years[16] and can now be predicted with high accuracy based on atomic structures of individual subunits.[17] Thus, identification of subcomplexes involved in assembly and disassembly processes has become technically feasible, then providing important clues to characterize intermediate states consisting of assembly and disassembly processes.

As for physicochemical analyses of mechanisms for transition between such intermediate states, atomistic simulations are practically useful at present; we could promptly find a feasible approach for the problem in a class of conventional free-energy profile (FEP) calculation protocols, e.g., umbrella sampling (US) combined with reaction pathway sampling techniques.[18,19] FEP such as potential of mean force (PMF) can provide quantitative insights into physicochemical mechanisms for biomacromolecule complex formation in terms of evaluated rate constants of configurational transition between two different intermediate states. Then, it could seem practically feasible to examine physicochemical mechanisms for multimeric biomacromolecule assembly and disassembly processes.

Nevertheless, there still remains one technical problem upon combining the subcomplex identification methods with FEP calculation protocols. They need a priori knowledge of a complete set of atomic coordinates for all subunits and subcomplexes at intermediate states in assembly or disassembly processes. It is not straightforward for molecular simulation approaches to provide the essential information. Even elementary events of assembly and disassembly processes (dimeric association and dissociation reactions, respectively) rarely occur within brute-force molecular dynamics (MD) simulations available of today.[20−24]

Reaction pathway sampling methods could be expected to overcome realistic restriction of computational resources although conventional ones are designed to work with a priori knowledge for a complete set of atomic coordinates of macrobiomolecules at each intermediate state.[25,26] It is still challenging for the above-mentioned state-of-the-art MS and SB approaches to provide atomic information sufficient for such pathway sampling simulations.

MS approaches identify subunits and subcomplex species as particles with molecular mass and net charge, so that their spatial resolution is bound to the macromolecular level. SB approaches directly predict a final docking pose and a pair of dissociated subcomplexes. It is therefore out of their methodological design to identify a configuration of a complete set of subunits, at each intermediate state involved in a pre-binding and post-unbound conditions.

Thus, a simulation scheme should be newly designed to solve the technical problem which prevents MD simulations from their application to assembly and disassembly processes associated with multiple biomacromolecules. We here propose a new reaction pathway sampling method. The theoretical framework we here employ is a hybrid Monte Carlo (hMC)/molecular dynamics (MD) scheme, since configuration sampling method in the MD phase is arbitrarily selected for specific research purposes, then resulting in flexible applications of this simulation framework to a variety of molecular phenomena in condensed phases.[27−30]

We can find research purpose-oriented flexibility within magnificent examples of hMC/MD simulation studies. Chemical bond formation and scission are empirically simulated by switching potential functions and using effective potential function in the frameworks of REDMOON[28] and KIMMDY,[30] respectively. Meanwhile, transform and relax sampling (TRS)[27] and hybrid neMD-MC[29] employ randomly assigned force bias and perturbed potential energy to accelerate transformation of biopolymer conformation, respectively.

In this method, we employed steered molecular dynamics (SMD) method to enhance nonequilibrium subprocesses composing assembly and disassembly processes, namely, inter-subunit association and dissociation reactions. This method has been widely used to simulate these reactions involved in dimeric biomacromolecule complexes.[19,31,32] Aiming additional enhancement of configuration sampling, we applied the configuration bias Monte Carlo (cbMC) scheme[33−35] to hybrid MC/MD simulation procedure, recalling the earlier hybrid MC/MD study resulting in several time enhancement of lipid bilayer configuration sampling[36] (see Sections SI-1 and SI-2 for further explanation for the design concept of our hcbMC/MD method).

To test our simulation method, referred to as hybrid configuration bias MC/MD (hcbMC/MD) method hereafter, we initially considered disassembly processes because of technical straightforwardness. The hcbMC/MD method is applicable to simulate disassembly processes with simply employing inter-subunit dissociation reactions without subunit–subunit docking pose prediction which is possibly needed to examine assembly processes associated with multiple subunits. In this study, the hcbMC/MD method was applied to simulate disassembly processes of serum amyloid P (SAP) component homo pentamer (Figure ), whose subcomplex species appearing in disassembly processes have been extensively studied in the previous MS studies.[13,37,38]

Figure 1

Molecular structures of serum amyloid P (SAP) component homo pentamer. (A) SAP pentamer with annotation for each subunit. (B) Two individual inter-subunit interaction surfaces of SAP subunit, highlighted by blue and red colors; S1 and S2 are depicted as representative of subunit pair. (C) Initial stages of SAP pentamer disassembly processes predicted by structural bioinformatics approach.[13] Reactions shown in (C) are supposed for any possible monomer and dimer.

Results and Discussion

Subcomplex Species and Disassembly Pathway Obtained from the hcbMC/MD Simulations

We performed five independent hcbMC/MD simulations for the SAP system. A hcbMC/MD cycle is repeated 500 times so that the total simulation time in the MD phase is 100 ns. The number of trial configuration M (see the Computational Methods section for the definition) was set to 10 as an attempt; choice of M could be further optimized to enhance sampling. Nonetheless, the sampling efficiency of our method is similar to that of the TRS method[27] and seems sufficient to promote disassembly processes.

A ring-shaped SAP pentamer is converted into a set of subcomplexes and subunits by the 500th cycle of hcbMC/MD simulations (Figure ). We find a subunit pair with nonspecific contacts, which were formed irrespective of inter-subunit interaction surfaces of SAP (see Figure B for comparison). Such nonspecific inter-subunit reassociation results from a periodic simulation box whose size is insufficient for apparent spatial separation among five SAP subunit. Then, we identified subcomplexes by excluding such a reassociated subunit pair that solely makes nonspecific inter-subunit contacts.

Figure 2

Configurations of SAP subunits at the last (500th) hcbMC/MD cycle, where each subunit pair has no native contacts (NC), thus being supposed to be dissociated. Two individual inter-subunit interaction surfaces of each subunit are highlighted by blue and red colors as in the case of Figure . The five individual simulations discussed in (A)–(E) are indexed by lowercase alphabets (a–e) hereafter.

Figure shows subcomplex species at each hcbMC/MD cycle for the homo-pentameric SAP system. In our hcbMC/MD simulations, any of subcomplex species (tetramer, trimer, and dimer) observed in the MS experiments[13] appears. MS experiments induce inter-subunit dissociation reactions using atomic collision with gas molecules, while hcbMC/MD simulations accelerate these reactions by employing repulsive forces acting between a subunit pair. Although our hcbMC/MD scheme employs dissociation mechanisms differing from those in the MS experimental method,[13] our simulations might finely evaluate relative SAP subunit interface strength similarly to the MS method.

Figure 3

Disassembly of serum amyloid P (SAP) component homo pentamer through hcbMC/MD cycles. The five individual simulations discussed in (A)–(E) correspond to those indexed by lowercase alphabets (a–e) in Figure . A distribution of SAP (sub)complex species is illustrated by color. The green and red arrows indicate the initial stage of SAP disassembly (generation of tetramer or trimer) and the first cycle where SAP pentamer undergoes a ring-opening event.

Besides, each of the five simulations results in decomposition into the five isolated SAP subunits. Meanwhile, the initial stage of SAP disassembly is classified into two distinct events: trimer–dimer formation and tetramer–monomer formation from the pentamer (Figure A,D,E and B,C, respectively; see Figure C). In Figure A, the trimer–dimer pair found at the 29th cycle returns to the pentamer form so that we suppose tetramer–monomer formation as the initial subcomplex disassembly. These two different disassembly events were similarly predicted by the structural bioinformatics approach in the earlier MS study.[13]

It is thus possible to say that the simulation results are consistent with those of the earlier study[13] with respect to SAP subcomplex species and the initial stage of SAP disassembly processes. Our reaction pathway sampling method could play a role complementary to MS and SB approaches in arrangement of complete sets of atomic coordinates of all subunits involved in a disassembly process.

Computation Performance of Hybrid Configuration Bias MC/MD Simulations

The disassembly of pentameric form of SAP can be efficiently simulated due to the methodological advantage of our hcbMC/MD approach. We can show this point using independent five unbiased NPT MD simulations (300 K, 1 bar) with simulation length of 200 ns in total, twice longer than those of hcbMC/MD simulations. As shown in Figure , the five SAP molecules retain each inter-subunit contacts during the 200 ns simulations, thus keeping to take a pentameric form in each of the unbiased MD simulations.

Figure 4

Averaged number of native contacts between SAP subunit pair. Average and statistical error values are calculated using 2000 snapshot structures obtained from each of the five unbiased 200 ns MD simulations (an error bar is buried in the square, then not being apparent in this panel). Each inter-subunit native contact is distinguished by different colors. Annotation of subunit pair (S1–S2, etc.) and indexes of individual simulation (a–e) are given as in the case of Figure and that of Figure , respectively. Statistical errors were estimated from standard error by considering 95% confidence interval.

Furthermore, each of the subunits retains folded configurations through hcbMC/MD cycles. Although we accelerated dissociation reactions using external forces, the values of root-mean-square deviation (RMSd) for SAP subunits are similar between hcbMC/MD simulations and unbiased 200 ns NPT MD simulations (Table ). This observation clearly shows that external forces applied to subunits in SMD simulations hardly affect their structures, thus additionally supporting applicability of our hcbMC/MD method to analyze physicochemical properties of proteins under physiological conditions.

Table 1

Root-Mean-Square Deviation (RMSd)a with Error Valueb for SAP Subunitc

	S1	S2	S3	S4	S5
hcbMC/MD	0.9 ± 0.2	0.9 ± 0.1	0.8 ± 0.1	0.8 ± 0.1	0.8 ± 0.1
unbiased MD	0.9 ± 0.2	0.8 ± 0.2	0.8 ± 0.2	0.8 ± 0.2	0.8 ± 0.3

RMSd in units of Å for each SAP subunit was averaged over each 500 hcbMC/MD cycles (100 ns in total) or 2000 unbiased MD cycle (200 ns in total), and averaged over the five simulations.

Error values denote CI95, derived from standard deviation over the five simulations.

SAP subunit is distinguished by combination of alphabetic letter and digit as in the case of Figure .

Free-Energy Profile of Ring-Opening Reaction of SAP Pentamer

Observing the hcbMC/MD trajectories, we can find a preliminary dissociation event in each of the five simulations. One of the five subunit interaction interfaces is broken, then resulting in ring-opened pentameric form of SAP (Figure ). A ring opening of the SAP pentamer occurs in a relatively early stage of hcbMC/MD cycles (Table ) and is followed by the initial steps of SAP pentamer disassembly process, that is, tetramer–monomer and trimer–dimer dissociation reactions (see Figure ).

Figure 5

Ring-opened SAP pentamer configuration obtained from hcbMC/MD simulation, indexed by “a” in Figure . The digits in (A) and (B) denote values of distance between centers of mass (COMs) of S3 and S4. The orange dotted line is depicted in the vicinity of broken subunit interaction interface. SAP subunit is distinguished by combination of alphabetic letter and digit as in the case of Figure .

Table 2

Geometrical Characterization of the SAP Pentamer at the Cycles Where the Initial Ring-Opening Reaction Occurs

simulation indexa	pair of SAP subunitb	cycle	CCSc (nm2)
a	S3–S4	20th	59.9
b	S1–S2	61st	59.8
c	S3–S4	11th	59.1
d	S2–S3	7th	59.3
e	S3–S4	10th	59.8

Each hcbMC/MD simulation is indexed by lowercase alphabet letters as in the case of Figure .

SAP subunit is distinguished by combination of alphabetic letter and digit as in the case of Figure .

CCS denotes collision cross section.

This ring-opening reaction has not been reported in the earlier mass spectroscopy studies,[13,37,38] although there are rooms to discuss the occurrence of the intermediate states under realistic experimental conditions. It is possible that MS observations have not detected such intermediate states due to their limitation for spatial resolution. Values of collision cross section (CCS) for these ring-opened SAP pentamer structures are around 59 nm2, being larger by ca. 3 nm2 than that calculated from the X-ray crystallographically resolved ring-forming pentamer.[39] A ring-opened SAP pentamer appearing in experimental conditions might be assigned into partially unfolded SAP pentamer in the context of mass spectroscopy.[37]

We then consider the possibility of experimental observation of this ring-opening reaction. High-speed time lapse atomic force microscopy (AFM) is now available to observe molecular dynamics of multimeric protein complexes with nanometer spatial resolution and hundreds of milliseconds temporal resolution.[40] This method could be employed to validate the presence of this reaction, if the reaction time scale is found within the observation period of AFM. Aiming to estimate the reaction time, we carried out free-energy profile calculations with a conventional protocol, umbrella sampling molecular dynamics (USMD) simulations combined with weighted histogram analyses method (WHAM).[41,42]

We here examined the initial ring-opening events appearing in the hcbMC/MD simulations due to the following observation for the trajectories we obtained. Ring-opened SAP pentamer does not necessarily lead to immediate progress to the initial dissociation steps but may transiently come back to ring-forming intact SAP pentamer (see Figures A and S3). A SAP structure fluctuating between the two pentameric forms might gradually lose inter-SAP subunit contacts; such partially broken contacts would give lower activation barrier of the ring-opening reaction so that additional free-energy calculations are probably needed to estimate a precedent free-energy loss.

We then arranged a set of initial atomic coordinates for USMD simulations using snapshot structures appearing at the cycle of initial ring-opening reaction and also those around the neighboring ones if needed. Each of these snapshot structures takes either ring-formed pentamer or ring-opened pentamer. A reaction coordinate is set to the distance between centers of mass (COMs) of subunits undergoing breakage of inter-subunit contacts (see Table for subunit pairs examined for PMF calculations and Section SI-3 in the Supporting Information for computational details of USMD simulations). The colored lines in Figure A–C show the PMFs derived from hcbMC/MD sampling-combined USMD simulations and the five independent hcbMC/MD simulations are distinguished by different colors.

Figure 6

Potential of mean force (PMF) of SAP pentamer ring-opening reaction (A)–(C), and reaction time evaluated with Eyring’s transition-state theory (D). Each hcbMC/MD simulation is indexed by lowercase alphabetical letter as in the case of Figure , and the corresponding PMF is distinguished from the remaining on panel by different colors (red: a; blue: b; green: c; purple: d; orange: e). In (A)–(C), PMFs calculated from the atomic coordinates obtained from the unbiased 10 ns NPT MD simulation, i.e., the initial structure for each of the hcbMC/MD simulations, are shown by gray lines. An error bar was calculated from standard deviation and denoted as CI95. SAP subunit pair is distinguished by combination of alphabetical letter and digit as in the case of Figure . In (A), the green line is found on the gray one. In (D), points above the red dotted borderline, netted with the orange square area, have the estimated reaction time falling within the observation period with AFM experiments.

We also calculated PMFs using the atomic coordinates obtained from the unbiased 10 ns NPT MD simulation, the initial structure for each of the hcbMC/MD simulations, as references for stable SAP pentamer configuration (see Section SI-4 in the Supporting Information for computational details of SMD-combined USMD simulations for the reference systems). These PMFs are shown by gray lines in Figure A–C.

These PMFs are shown by gray lines in Figure A–C. Each of PMFs has one activation barrier and appears to be up-hill (Figure A–C). A steep change of the slope of PMF is found between 37 and 40 Å on the reaction coordinates (see also Figure A,B for configurational change of SAP pentamer). Since we here characterize the reaction mechanism using one reaction coordinate, considerable deviations among the values of activation barrier probably indicate the presence and the effect of other degrees of freedom orthogonal to the distance between COMs of subunits, which have substantial influence on the reaction process.

The significant decrease of activation barrier found in the simulation indexed by “b” (referred to as “Sim-b” hereafter) could be explained by difference in the number of atomic contacts between a subunit pair which retains contacts of inter-subunit interaction surfaces. Figure shows the number of atomic contacts between each subunit pair for the five hcbMC/MD simulations, where the snapshot structures are obtained from the simulation cycle at which we observed the initial ring-opening reaction. With regard to inter-subunit atomic contacts retained in ring-opened SAP pentamer, we can find a smaller number of atomic contacts in S5–S1 interface in Sim-b, that is, 8. Meanwhile, the remaining three are 15 or larger. It is of note that S5 is neighboring to S1, which undergoes ring-opening reaction in Sim-b (Table ).

Figure 7

Atomic contacts between subunit pair for snapshot structure at hcbMC/MD cycle for the initial ring-opening reaction. Annotation of SAP subunit pair (S1–S2, etc.) and indexes of individual simulation (a–e) are given as in the case of Figure and that of Figure , respectively.

From the viewpoint of reaction energetics, stepwise breakage of subunit interface seems to more easily occur rather than the simultaneous break of multiple subunit interfaces. Recalling the earlier simulation study on stepwise dissociation reaction of protein dimer,[22] we can assume that this ring-opening reaction proceeds in a similar manner; allosteric losing of inter-subunit contacts is followed by dissociation between a subunit pair.

We suppose that breakage of salt bridge formed between Lys117 and Asp42 in S5–S1 interface (Figure A–C) is particularly important to weaken configurational restraints between S1 and S5 in the ring-opened SAP pentamer in Sim-b. This salt bridge breakage seems to allow libration motion of S1 on S5 and enhances dissociation of S1 from S2. The effects of Lys117 and Asp42 on SAP pentamer formation have not been reported to our best knowledge. While considering significant decrease of the activation barrier for Sim-b, mutation on either or both of these two charged residues, or strongly acidic and basic conditions, might repress SAP oligomerization. SAP was identified as pathological deposits in 1967 and has been widely studied in the field of amyloid- and immune-related diseases.[43,44] These two residues might make a theraphetic target in aim of regulation of SAP pentamer formation.

Figure 8

Positional relation between Asp42 and Lys117 at the subunit interface. (A) The initial structure of hcbMC/MD simulations. (B) Structure at the timing of initial ring opening in the simulation indexed with (B). (C) Distances between Cγ in Asp42 and Nζ in Lys117 for each interface. SAP subunit is distinguished by combination of alphabetic letter and digit as in the case of Figure . In (A), the blue dotted lines denote hydrogen bonding between Lys117 and Asp42. In panel (C), the red dotted line shows the level of 4 Å. Each hcbMC/MD simulation is indexed by lowercase alphabet letters (a–e) as in the case of Figure , and the index “i” denotes the atomic coordinates obtained from the unbiased 10 ns NPT MD simulation, i.e., the initial structure for each of the hcbMC/MD simulations.

According to the above discussion, we can speculate that atomistic events occurring in the hcbMC/MD simulation, such as preliminary deformation at the vicinal SAP subunit interface, has influence on activation barrier of the ring-opening reaction. This speculation is supported from comparison between the two free-energy profiles; one for ring-opening simulated from an hcbMC/MD trajectory, other for ring-opening simulated with the initial structure of hcbMC/MD, that is, reference to stable SAP pentamer configuration. We can find substantial free-energy decrease of 11 kcal/mol in Sim-b and that of 6 kcal/mol in the simulation indexed by “d” at 50 Å on each reaction coordinate (see Figure B,C, respectively).

In particular, compared with the other four simulations, the change for Sim-b is so large as to make the time scale of the ring-opening reaction 107-fold smaller than those obtained from the other four simulations (see Figure D). The above observation indicates physicochemical importance of the preliminary deformation. The ring-opening event observed in Sim-b might be essentially a two-step reaction accompanying the preliminary deformation. Besides, the preliminary deformation itself is possibly a time-consuming process, which is characterized by free-energy profile with barrier of a substantial height, thus being the rate-determining step of ring-opening reaction pathway observed in Sim-b.

We have here discussed the reaction mechanism in terms of one reaction coordinate, distance between COMs of a subunit pair involved in interface breakage. The considerable deviations among the one-dimensional PMFs are probably due to a degree of freedom orthogonal to the distance between COMs of subunits, which is associated with preliminary deformation at the vicinal SAP subunit interface. Then, we could refine physicochemical characterization of the reaction mechanism by employing an additional reaction coordinate, e.g., the number of allosteric inter-subunit contacts at the vicinal interface. The free-energy profiles discussed in Figure A–C might be mapped as individual pathways on a multidimensional free-energy landscape. Similarly, the time scale of preliminary deformation could also be elucidated quantitatively and the actual reaction time of ring-opening pathway in Sim-b could be obtained. Nonetheless, further effort to address these problems is beyond the scope of this study, thus being left for future studies.

While considering physicochemical properties of the free-energy profiles, it could be supposed that the ring-opening reaction observed in our simulation is experimentally testable. Among the five hcbMC/MD-derived PMFs, the simulation indexed by “c” shows the highest activation barrier of 23.7 kcal/mol (see the blue line in Figure A), and the corresponding reaction time is 8.2 h according to eq . This estimated reaction time may fall within observation periods of time lapse AFM, hundreds of milliseconds or longer, so that occurrence of this ring-opening reaction could be validated by employing AFM experiments. The above discussion leads to such a speculation that a part of SAP pentamer might have ring-opening forms during the AMF observation period.

It is noted that the above speculation is similarly obtained even if we use Kramers’ transition-state theory. Kramers’ theory possibly give a longer reaction time scale than Eyring’s one,[45] although this does not essentially change the speculation; ring-opening reactions are still supposed to occur within the AFM observation time scale.

Recalling observation period of time lapse AFM experiments, we then classify the five ring-opening events into two groups; that for Sim-b and those for the other four. The reaction time obtained from Sim-b is 39 ns, significantly smaller than those obtained from the other four simulations (Figure B). The reaction time of 39 ns means that a ring-opening reaction can spontaneously proceed under thermal fluctuation along this reaction pathway. Such faster molecular dynamics would be monitored using fluorescence resonance energy transfer (FRET), which has been extensively used to analyze dynamics of protein complexes (additional remarks on the usage of FRET are given in Section SI-5 in the Supporting Information).[46,47]

Concluding Remarks

We have here proposed a new reaction pathway sampling method, hcbMC/MD, aiming to simulate complicated assembly and disassembly processes involving multimeric biomacromolecules within practically accessible computation time.

Our method was applied to disassembly simulations of serum amyloid P (SAP) component homo pentamer. Subcomplex species and disassembly pathway we obtained from hcbMC/MD simulations are consistent with the earlier observations. Our simulation approach could play a role complementary to mass spectrometry and structural bioinformatics approaches with regard to arrangement of complete sets of atomic coordinates of a molecular system undergoing disassembly processes.

Regarding observations obtained from mass spectroscopy and structural bioinformatics approaches as molecular and structural landmarks, our atomistic simulation approaches could provide deeper insights into physicochemical mechanism for biomacromolecular assembly and disassembly processes. In fact, we observed a novel dissociation event, ring-opening reaction of SAP pentamer. Employing free-energy calculation combined with our hcbMC/MD trajectories, we moreover obtained experimentally testable observations on (1) reaction time of the ring-opening reaction and (2) importance of Asp42 and Lys117 for stable formation of SAP oligomer.

Considering physicochemical complexity of atomistic multicomponent systems and practical restriction of computational resources, a configuration ensemble generated by any enhanced sampling method is probably a subset of an exact configurational space. Although this viewpoint often raises such a concern that enhanced sampling simulations for biomacromolecules pick up biologically insignificant configurations, we suppose that our hcbMC/MD method is practically well designed to sample various kinds of macromolecular configurations associated with multimeric biomacromolecule disassembly pathway (metastable complex species, subcomplex species, and isolated subunits) and examine the mechanisms at the atomic level (an additional remark on our hcbMC/MD-derived configuration ensemble is given in Section SI-6 in the Supporting Information).

As shown above, our pathway sampling simulations succeeded to convert SAP pentamer into the five SAP subunits within realistic computational time. SAP subcomplex species appearing through the disassembly processes were consistent with those deduced from the MS experiments. Furthermore, we obtained experimentally testable predictions for the disassembly processes and mechanisms via free-energy calculations. We then suppose that our hcbMC/MD approach opens a new avenue to study physicochemical mechanisms of complicated biological processes associated with multimeric biomacromolecules in the cell.

Computational Methods

Setup of Serum Amyloid P Component System

We used atomic coordinates of serum amyloid P (SAP) component homo pentamer (Protein Data Bank (PDB) entry: 4AVS).[39] Two calcium ions bound to SAP subunits were retained in the atomic coordinates of this SAP system. Nε protonation state was employed for each of histidine residues, and all carboxyl groups in aspartate and glutamate residues were set to the deprotonated state, where we considered physiological pH condition in the cell. The SAP structure was solvated in the rectangular box with 49 305 water molecules and was electrically neutralized by adding 10 Cl– ion molecules.

To calculate the forces acting among atoms, AMBER force field 14SB,[48] TIP3P water model,[49,50] and JC ion parameters adjusted for the TIP3P water model[51,52] were employed for amino acid residues, water molecules, and ions, respectively. In addition, the force field parameter sets developed by Bradbrook was applied for divalent calcium ion.[53] Molecular modeling of the SAP system was performed using the LEaP modules in AmberTools 17 package.[54] Molecular mechanics (MM) and molecular dynamics (MD) simulations were performed under the periodic boundary condition with graphics processing unit (GPU)-version PMEMD module in AMBER 17 package[54] based on SPFP algorithm[55] with NVIDIA GeForce GTX1080 Ti. The solvated structure was relaxed through 10 ns NPT MD simulations, then being employed as the initial atomic coordinates for the following hybrid configuration bias MC/MD simulations. Further computational details are described in Sections SI-7 and SI-8 in the Supporting Information.

Hybrid Configuration Bias MC/MD Simulation Scheme

The simulation scheme is designed based on a general configuration bias MC method[33] as follows. We here employed unbiased MD and SMD simulations to generate old and new sets of atomic configurations of molecular systems, respectively. Inter-subunit dissociation reactions are mainly accelerated with SMD simulations. Detailed simulation conditions of MD and SMD-based configuration generation are given in the next subsection.

Generate M trial configurations {a1, a2, ..., a} by an unbiased MD simulation and calculate bias energy ubias, which is used for biasing configuration selection and will be specified later, for each configuration.

Assign the last snapshot obtained from the simulation (a) as an old configuration, denoted by X and define the Rosenbluth factor[33]β denotes inverse of kBT, where kB and T are the Boltzmann constant and system temperature, respectively. In the following step 3, new configurations are generated by starting from a so that we assigned this configuration to x by recalling conventional cbMC schemes.[33]

Generate M trial configurations {b1, b2, ..., b} by a steered MD simulation starting from the old configuration x and calculate bias energy ubias for each configuration.

Define the Rosenbluth factorand select one configuration among {b1, b2, ..., b}, denoted by x, with a probability

The configuration change is accepted with a probabilitywhere x and U(x) denote a configuration and the native potential energy function of the system with configuration of x, respectively.

Considering that configurations of the system are generated by moving all atoms simultaneously with the molecular dynamics method, we used the Rosenbluth factor in the forms given in eqs and 2. Monte Carlo trial with this type of configuration sampling is particularly referred to as orientation bias Monte Carlo.[33]

As for the right-hand side of eq , we assume that the probability of generating forward move (x → x) is equal to that of generating backward move (x ← x). This assumption is employed here for simplification of simulation implementation but not trivial in general. Our hcbMC/MD method could work for sampling various kinds of macromolecular configurations, including metastable complex species, subcomplex species, and isolated subunits, by employing the SMD method, while it is uncertain to exactly construct a canonical ensemble defined by an unperturbed Hamiltonian. Then, with the aim of validating physicochemical observations obtained from an hcbMC/MD trajectory, it is a practical option to perform additional umbrella sampling (US) MD simulations with a set of snapshot structures involved in an inter-subunit dissociation event during the trajectory.

A multimeric protein complex disassembly process should be accompanied by gradual breakage of inter-subunit native contacts (NC). Then, we defined the bias energy as a function of the total number of inter-subunit NC formed in a multimeric protein complex (nNC), where the initial atomic coordinates are used as the reference structure:

The coefficient k is set to β–1 at 300 K, that is, 0.6 kcal/mol by assuming energetic stabilization obtained from hydrophobic atomic contacts.[56] A set of native contacts at subunit binding interface was defined using the initial atomic coordinates of SAP pentamer for hcbMC/MD simulations. All hcbMC/MD simulations were performed by employing an in-house MC/MD simulation interface.

Configuration Generation at Each hcbMC/MD Cycle Using Unbiased and Steered MD Simulations

Configurations of SAP system, denoted by {a1, a2, ..., a}, were sampled using an unbiased NPT MD simulation (300 K, 1 bar). Meanwhile, configurations of the system undergoing a partial inter-subunit dissociation, denoted by {b1, b2, ..., b}, were sampled using a steered molecular dynamics (SMD) simulation under the NPT condition (300 K, 1 bar), starting from the snapshot structure obtained from the unbiased MD simulation. The simulation length is 100 ps for each of MD and SMD simulations. A set of atomic coordinates was recorded by every 10 ps and was indexed in ascending order to prepare {a1, a2, ..., a} and {b1, b2, ..., b}; M is set to 10 in this simulation condition. The system temperature and pressure were regulated by Langevin thermostat with 1 ps–1 collision coefficient, and Monte Carlo barostat with the attempt of system volume change by every 100 steps, respectively.

In every hcbMC/MD cycle, we randomly chose a pair of monomers in the 10 (=5C2, derived from a combination of two subunits among the five) candidates. We employed the distance between COMs of the chosen subunit pair as the reaction coordinate of SMD simulation. A COM for an SAP subunit was calculated using all of the 204 Cα atoms.

In an SMD simulation, a target value of the distance was set to d0 + Δd, where d0 and Δd are the initial value of the distance at the cycle and a random integer in the range of 8–12, respectively. The harmonic potential with the force constant of 10 kcal/(mol Å2) was imposed on the reaction coordinate. In each cycle, a reaction coordinate for inter-SAP subunit dissociation reaction and a random seed for Langevin thermostat were changed randomly. Following the above simulation conditions, an unbiased MD simulation was also designed to evaluate the sampling performance of hcbMC/MD simulation, where the 100 ps MD simulation cycle was repeated 2000 times. Additional computational details and the corresponding unbiased MD simulation procedure are described in Sections SI-7 and SI-9 in the Supporting Information, respectively.

Umbrella Sampling Molecular Dynamics Simulations

Certain numbers of snapshot structures were extracted from a hcbMC/MD trajectory and employed for USMD windows. Following temperature relaxation simulations, 5 ns NVT USMD simulations (300 K) were performed for each of the USMD windows (see Table S1 in the Supporting Information for details). The system temperature was regulated using Langevin thermostat with 1 ps–1 collision coefficient. Each of the last 3 ns USMD trajectories was used to construct a PMF (see Figures S1 and S2 for convergence of PMF curves). A reaction coordinate is the distance between centers of mass (COM) of the subunit pair undergoing dissociation reaction, since an inter-subunit COM distance has been widely used as reaction coordinate to characterize association and dissociation reactions.[19,32] A set of USMD simulations was repeated eight times. Besides, we calculated PMFs of stable SAP pentamer by employing the structure obtained from the unbiased 10 ns NPT MD simulation, as the initial atomic coordinates for each of the hcbMC/MD simulations. Further technical details are given in the Supporting Information (see Sections SI-3 and SI-4).

Trajectory Analyses

Root-mean-square deviation (RMSd), and native and non-native atomic contacts were analyzed with the cpptraj module in AmberTools 17 package.[54] We calculated RMSd and native contacts using the Cα atoms and the nonhydrogen atoms in the initial atomic coordinates for the hcbMC/MD simulations, respectively. The distance criterion for atomic contacts was set to 3.5 Å. Inter-subunit dissociation and reassociation conditions are defined as follows: Nnative = 0 and Nnon-native = 0; Nnon-native ≥ 0 in Nnative > 0, respectively. Nnative and Nnon-native denote the number of native atomic contacts and that of non-native atomic contacts made between a pair of SAP subunit, respectively. Collision cross section (CCS) of SAP protein was analyzed using IMPACT program.[57]

PMF was calculated with weighted histogram analysis method (WHAM)[41,42] using each set of USMD trajectories. Statistical errors of PMF values, σPMF(ξ), were estimated by bootstrapped sampling:[58]Here, Nb, ξ, and Wb,(ξ) denote the number of bootstrapped sampling, the reaction coordinate, and the value of kth bootstrapped PMF at each point of ξ, respectively. ⟨Wb(ξ)⟩ is the average over all ⟨Wb,(ξ)⟩, where the value of Nb is set to 200.

Reaction rate, kTST, is estimated using Eyring’s transition-state theory:where ΔF†, h, kB, and T denote the activation barrier height, Planck’s constant, Boltzmann constant, and the temperature of the system, respectively. The estimation by employing formula is supposed to be an upper bound of the reaction rate.[45]

Reaction time scale, τTST, is defined as the inverse of kTST. ΔF† is defined as F(ξ0′) – F(ξ0), where the value of PMF has a local minimum at ξ0, and a maximum at ξ0′, which is greater than ξ0.

Molecular structures were illustrated using Visual Molecular Dynamics (VMD).[59] Error bars were calculated from standard error that indicate 95% confidence interval.