Structure of the Dnmt1 Reader Module Complexed with a Unique Two-Mono-Ubiquitin Mark on Histone H3 Reveals the Basis for DNA Methylation Maintenance.

The proper location and timing of Dnmt1 activation are essential for DNA methylation maintenance. We demonstrate here that Dnmt1 utilizes two-mono-ubiquitylated histone H3 as a unique ubiquitin mark for its recruitment to and activation at DNA methylation sites. The crystal structure of the replication foci targeting sequence (RFTS) of Dnmt1 in complex with H3-K18Ub/23Ub reveals striking differences to the known ubiquitin-recognition structures. The two ubiquitins are simultaneously bound to the RFTS with a combination of canonical hydrophobic and atypical hydrophilic interactions. The C-lobe of RFTS, together with the K23Ub surface, also recognizes the N-terminal tail of H3. The binding of H3-K18Ub/23Ub results in spatial rearrangement of two lobes in the RFTS, suggesting the opening of its active site. Actually, incubation of Dnmt1 with H3-K18Ub/23Ub increases its catalytic activity in vitro. Our results therefore shed light on the essential role of a unique ubiquitin-binding module in DNA methylation maintenance.

Structures and mechanisms of actin ATP hydrolysis.

The major cytoskeleton protein actin undergoes cyclic transitions between the monomeric G-form and the filamentous F-form, which drive organelle transport and cell motility. This mechanical work is driven by the ATPase activity at the catalytic site in the F-form. For deeper understanding of the actin cellular functions, the reaction mechanism must be elucidated. Here, we show that a single actin molecule is trapped in the F-form by fragmin domain-1 binding and present their crystal structures in the ATP analog-, ADP-Pi-, and ADP-bound forms, at 1.15-Å resolutions. The G-to-F conformational transition shifts the side chains of Gln137 and His161, which relocate four water molecules including W1 (attacking water) and W2 (helping water) to facilitate the hydrolysis. By applying quantum mechanics/molecular mechanics calculations to the structures, we have revealed a consistent and comprehensive reaction path of ATP hydrolysis by the F-form actin. The reaction path consists of four steps: 1) W1 and W2 rotations; 2) PG-O3B bond cleavage; 3) four concomitant events: W1-PO3- formation, OH- and proton cleavage, nucleophilic attack by the OH- against PG, and the abstracted proton transfer; and 4) proton relocation that stabilizes the ADP-Pi-bound F-form actin. The mechanism explains the slow rate of ATP hydrolysis by actin and the irreversibility of the hydrolysis reaction. While the catalytic strategy of actin ATP hydrolysis is essentially the same as those of motor proteins like myosin, the process after the hydrolysis is distinct and discussed in terms of Pi release, F-form destabilization, and global conformational changes.

Conserved gatekeeper methionine regulates the binding and access of kinase inhibitors to ATP sites of MAP2K1, 4, and 7: Clues for developing selective inhibitors.

Mitogen-activated protein kinase kinases (MAP2Ks) 1, 4, and 7 are potential targets for treating various diseases. Here, we solved the crystal structures of MAP2K1 and MAP2K4 complexed with covalent inhibitor 5Z-7-oxozeaenol (5Z7O). The elucidated structures showed that 5Z7O was non-covalently bound to the ATP binding site of MAP2K4, while it covalently attached to cysteine at the DFG-1 position of the deep ATP site of MAP2K1. In contrast, we previously showed that 5Z7O covalently binds to MAP2K7 via another cysteine on the solvent-accessible edge of the ATP site. Structural analyses and molecular dynamics calculations indicated that the configuration and mobility of conserved gatekeeper methionine located at the central ATP site regulated the binding and access of 5Z7O to the ATP site of MAP2Ks. These structural features provide clues for developing highly potent and selective inhibitors against MAP2Ks. Abbreviations: ATP, adenosine triphosphate; FDA, Food and Drug Administration; MAP2Ks, mitogen-activated protein kinase kinases; MD, molecular dynamics; NSCLC, non-small cell lung cancer; 5Z7O, 5Z-7-oxozeaenol; PDB, protein data bank; RMSD, root-mean-square deviation.

Binding and Unbinding Pathways of Peptide Substrates on the SARS-CoV-2 3CL Protease.

Based on many crystal structures of ligand complexes, much study has been devoted to understanding the molecular recognition of SARS-CoV-2 3C-like protease (3CLpro), a potent drug target for COVID-19. In this research, to extend this present static view, we examined the kinetic process of binding/unbinding of an eight-residue substrate peptide to/from 3CLpro by evaluating the path ensemble with the weighted ensemble simulation. The path ensemble showed the mechanism of how a highly flexible peptide folded into the bound form. At the early stage, the dominant motion was the diffusion on the protein surface showing a broad distribution, whose center was led into the cleft of the chymotrypsin fold. We observed a definite sequential formation of the hydrogen bonds at the later stage occurring in the cleft, initiated between Glu166 (3CLpro) and P3_Val (peptide), followed by binding to the oxyanion hole and completed by the sequence-specific recognition at P1_Gln.

A cutoff-based method with charge-distribution-data driven pair potentials for efficiently estimating electrostatic interactions in molecular systems.

We introduce a simple cutoff-based method for precise electrostatic energy calculations in the molecular dynamics (MD) simulations of point-particle systems. Our method employs a theoretically derived smooth pair potential function to define electrostatic energy, offering stability and computational efficiency in MD simulations. Instead of imposing specific physical conditions, such as dielectric environments or charge neutrality, we focus on the relationship represented by a single summation formula of charge-weighted pair potentials. This approach allows an accurate energy approximation for each particle, enabling a straightforward error analysis. The resulting particle-dependent pair potential captures the charge distribution information, making it suitable for heterogeneous systems and ensuring an enhanced accuracy through distant information inclusion. Numerical investigations of the Madelung constants of crystalline systems validate the method's accuracy.

Flexibility and Cell Permeability of Cyclic Ras-Inhibitor Peptides Revealed by the Coupled Nosé-Hoover Equation.

Quantifying the cell permeability of cyclic peptides is crucial for their rational drug design. However, the reasons remain unclear why a minor chemical modification, such as the difference between Ras inhibitors cyclorasin 9A5 and 9A54, can substantially change a peptide's permeability. To address this question, we performed enhanced sampling simulations of these two 11-mer peptides using the coupled Nosé-Hoover equation (cNH) we recently developed. The present cNH simulations realized temperature fluctuations over a wide range (240-600 K) in a dynamic manner, allowing structural samplings that were well validated by nuclear Overhauser effect measurements. The derived structural ensembles were comprehensively analyzed by all-atom structural clustering, mapping the derived clusters onto principal components (PCs) that characterize the cyclic structure, and calculating cluster-dependent geometric and chemical properties. The planar-open conformation was dominant in aqueous solvent, owing to inclusion of the Trp side chain in the main-chain ring, while the compact-closed conformation, which favors cell permeation due to its compactness and high polarity, was also accessible. Conformation-dependent cell permeability was observed in one of the derived PCs, demonstrating that decreased cell permeability in 9A54 is due to the high free energy barrier separating the two conformations. The origin of the change in free energy surface was determined to be loss of flexibility in the modified residues 2-3, resulting from the increased bulkiness of their side chains. The derived molecular mechanism of cell permeability highlights the significance of complete structural dynamics surveys for accelerating drug development with cyclic peptides.

Conformational change of a biomolecule studied by the weighted ensemble method: Use of the diffusion map method to extract reaction coordinates.

We simulate the nonequilibrium ensemble dynamics of a biomolecule using the weighted ensemble method, which was introduced in molecular dynamics simulations by Huber and Kim and further developed by Zuckerman and co-workers. As the order parameters to characterize its conformational change, we here use the coordinates derived from the diffusion map (DM) method, one of the manifold learning techniques. As a concrete example, we study the kinetic properties of a small peptide, chignolin in explicit water, and calculate the conformational change between the folded and misfolded states in a nonequilibrium way. We find that the transition time scales thus obtained are comparable to those using previously employed hydrogen-bond distances as the order parameters. Since the DM method only uses the 3D Cartesian coordinates of a peptide, this shows that the DM method can extract the important distance information of the peptide without relying on chemical intuition. The time scales are compared well with the previous results using different techniques, non-Markovian analysis and core-set milestoning for a single long trajectory. We also find that the most significant DM coordinate turns out to extract a dihedral angle of glycine, and the previously studied relaxation modes are well correlated with the most significant DM coordinates.

Multiscale enhanced sampling of glucokinase: Regulation of the enzymatic reaction via a large scale domain motion.

Enhanced sampling yields a comprehensive structural ensemble or a free energy landscape, which is beyond the capability of a conventional molecular dynamics simulation. Our recently developed multiscale enhanced sampling (MSES) method employs a coarse-grained model coupled with the target physical system for the efficient acceleration of the dynamics. MSES has demonstrated applicability to large protein systems in solution, such as intrinsically disordered proteins and protein-protein and protein-ligand interactions. Here, we applied the MSES simulation to an important drug discovery target, glucokinase (GCK), to elucidate the structural basis of the positive cooperativity of the enzymatic reaction at an atomistic resolution. MSES enabled us to compare two sets of the free energy landscapes of GCK, for the glucose-bound and glucose-unbound forms, and thus demonstrated the drastic change of the free energy surface depending on the glucose concentration. In the glucose-bound form, we found two distinct basins separated by a high energy barrier originating from the domain motion and the folding/unfolding of the α13 helix. By contrast, in the glucose-unbound form, a single flat basin extended to the open and super-open states. These features illustrated the two distinct phases achieving the cooperativity, the fast reaction cycle staying in the closed state at a high glucose concentration and the slow cycle primarily in the open/super-open state at a low concentration. The weighted ensemble simulations revealed the kinetics of the structural changes in GCK with the synergetic use of the MSES results; the rate constant of the transition between the closed state and the open/super-open states, kC/O = 1.1 ms-1, is on the same order as the experimental catalytic rate, kcat = 0.22 ms-1. Finally, we discuss the pharmacological activities of GCK activators (small molecular drugs modulating the GCK activity) in terms of the slight changes in the domain motion, depending on their chemical structures as regulators. The present study demonstrated the capability of the enhanced sampling and the associated kinetic calculations for understanding the atomistic structural dynamics of protein systems in physiological environments.

Exploring Configuration Space and Path Space of Biomolecules Using Enhanced Sampling Techniques-Searching for Mechanism and Kinetics of Biomolecular Functions.

To understand functions of biomolecules such as proteins, not only structures but their conformational change and kinetics need to be characterized, but its atomistic details are hard to obtain both experimentally and computationally. Here, we review our recent computational studies using novel enhanced sampling techniques for conformational sampling of biomolecules and calculations of their kinetics. For efficiently characterizing the free energy landscape of a biomolecule, we introduce the multiscale enhanced sampling method, which uses a combined system of atomistic and coarse-grained models. Based on the idea of Hamiltonian replica exchange, we can recover the statistical properties of the atomistic model without any biases. We next introduce the string method as a path search method to calculate the minimum free energy pathways along a multidimensional curve in high dimensional space. Finally we introduce novel methods to calculate kinetics of biomolecules based on the ideas of path sampling: one is the Onsager⁻Machlup action method, and the other is the weighted ensemble method. Some applications of the above methods to biomolecular systems are also discussed and illustrated.

Energy landscape of all-atom protein-protein interactions revealed by multiscale enhanced sampling.

Protein-protein interactions are regulated by a subtle balance of complicated atomic interactions and solvation at the interface. To understand such an elusive phenomenon, it is necessary to thoroughly survey the large configurational space from the stable complex structure to the dissociated states using the all-atom model in explicit solvent and to delineate the energy landscape of protein-protein interactions. In this study, we carried out a multiscale enhanced sampling (MSES) simulation of the formation of a barnase-barstar complex, which is a protein complex characterized by an extraordinary tight and fast binding, to determine the energy landscape of atomistic protein-protein interactions. The MSES adopts a multicopy and multiscale scheme to enable for the enhanced sampling of the all-atom model of large proteins including explicit solvent. During the 100-ns MSES simulation of the barnase-barstar system, we observed the association-dissociation processes of the atomistic protein complex in solution several times, which contained not only the native complex structure but also fully non-native configurations. The sampled distributions suggest that a large variety of non-native states went downhill to the stable complex structure, like a fast folding on a funnel-like potential. This funnel landscape is attributed to dominant configurations in the early stage of the association process characterized by near-native orientations, which will accelerate the native inter-molecular interactions. These configurations are guided mostly by the shape complementarity between barnase and barstar, and lead to the fast formation of the final complex structure along the downhill energy landscape.

Minimum free energy path of ligand-induced transition in adenylate kinase.

Large-scale conformational changes in proteins involve barrier-crossing transitions on the complex free energy surfaces of high-dimensional space. Such rare events cannot be efficiently captured by conventional molecular dynamics simulations. Here we show that, by combining the on-the-fly string method and the multi-state Bennett acceptance ratio (MBAR) method, the free energy profile of a conformational transition pathway in Escherichia coli adenylate kinase can be characterized in a high-dimensional space. The minimum free energy paths of the conformational transitions in adenylate kinase were explored by the on-the-fly string method in 20-dimensional space spanned by the 20 largest-amplitude principal modes, and the free energy and various kinds of average physical quantities along the pathways were successfully evaluated by the MBAR method. The influence of ligand binding on the pathways was characterized in terms of rigid-body motions of the lid-shaped ATP-binding domain (LID) and the AMP-binding (AMPbd) domains. It was found that the LID domain was able to partially close without the ligand, while the closure of the AMPbd domain required the ligand binding. The transition state ensemble of the ligand bound form was identified as those structures characterized by highly specific binding of the ligand to the AMPbd domain, and was validated by unrestrained MD simulations. It was also found that complete closure of the LID domain required the dehydration of solvents around the P-loop. These findings suggest that the interplay of the two different types of domain motion is an essential feature in the conformational transition of the enzyme.

Multiscale enhanced path sampling based on the Onsager-Machlup action: application to a model polymer.

We propose a novel path sampling method based on the Onsager-Machlup (OM) action by generalizing the multiscale enhanced sampling technique suggested by Moritsugu and co-workers [J. Chem. Phys. 133, 224105 (2010)]. The basic idea of this method is that the system we want to study (for example, some molecular system described by molecular mechanics) is coupled to a coarse-grained (CG) system, which can move more quickly and can be computed more efficiently than the original system. We simulate this combined system (original + CG system) using Langevin dynamics where different heat baths are coupled to the two systems. When the coupling is strong enough, the original system is guided by the CG system, and is able to sample the configuration and path space with more efficiency. We need to correct the bias caused by the coupling, however, by employing the Hamiltonian replica exchange, where we prepare many path replicas with different coupling strengths. As a result, an unbiased path ensemble for the original system can be found in the weakest coupling path ensemble. This strategy is easily implemented because a weight for a path calculated by the OM action is formally the same as the Boltzmann weight if we properly define the path "Hamiltonian." We apply this method to a model polymer with Asakura-Oosawa interaction, and compare the results with the conventional transition path sampling method.

Disorder-to-order transition of an intrinsically disordered region of sortase revealed by multiscale enhanced sampling.

Molecular functions of intrinsically disordered proteins (IDPs) or intrinsically disordered regions (IDRs), such as molecular recognition and cellular signaling, are ascribed to dynamic changes in the conformational space in response to binding of target molecules. Sortase, a transpeptitase in Gram-positive bacteria, has an IDR in a loop which undergoes a disordered-to-ordered transition (called "disordered loop"), accompanying a tilt of another loop ("dynamic loop"), upon binding of a signal peptide and a calcium ion. In this study, all-atom conformational ensembles of sortase were calculated for the four different binding states (with/without the peptide and with/without a calcium ion) by the multiscale enhanced sampling (MSES) simulation to examine how the binding of the peptide and/or calcium influences the conformational ensemble. The MSES is a multiscale and multicopy simulation method that allows an enhanced sampling of the all-atom model of large proteins including explicit solvent. A 100 ns MSES simulation of the ligand-free sortase using 20 replicas (in total 2 µs) demonstrated large flexibility in both the disordered and dynamic loops; however, their distributions were not random but had a clear preference which populates the N-terminal part of the disordered loop near the bound form. The MSES simulations of the three binding states clarified the allosteric mechanism of sortase: the N- and C-terminal parts of the disordered loop undergo a disorder-to-order transition independently of each other upon binding of the peptide and a calcium ion, respectively; however, upon binding of both ligands, the two parts work cooperatively to stabilize the bound peptide.

REACH coarse-grained simulation of a cellulose fiber.

A molecular level understanding of the structure, dynamics and mechanics of cellulose fibers can aid in understanding the recalcitrance of biomass to hydrolysis in cellulosic biofuel production. Here, a residue-scale REACH (Realistic Extension Algorithm via Covariance Hessian) coarse-grained force field was derived from all-atom molecular dynamics (MD) simulations of the crystalline Iß cellulose fibril. REACH maps the atomistic covariance matrix onto coarse-grained elastic force constants. The REACH force field was found to reproduce the positional fluctuations and low-frequency vibrational spectra from the all-atom model, allowing elastic properties of the cellulose fibril to be characterized using the coarse-grained force field with a speedup of >20 relative to atomistic MD on systems of the same size. The calculated longitudinal/transversal Young's modulus and the velocity of sound are in agreement with experiment. The persistence length of a 36-chain cellulose microcrystal was estimated to be ~380 µm. Finally, the normal-mode analysis with the REACH force field suggests that intrinsic dynamics might facilitate the deconstruction of the cellulose fibril from the hydrophobic surface.

Functional dynamics of SARS-CoV-2 3C-like protease as a member of clan PA.

SARS-CoV-2 3C-like protease (3CLpro), a potential therapeutic target for COVID-19, consists of a chymotrypsin fold and a C-terminal α-helical domain (domain III), the latter of which mediates dimerization required for catalytic activation. To gain further understanding of the functional dynamics of SARS-CoV-2 3CLpro, this review extends the scope to the comparative study of many crystal structures of proteases having the chymotrypsin fold (clan PA of the MEROPS database). First, the close correspondence between the zymogen-enzyme transformation in chymotrypsin and the allosteric dimerization activation in SARS-CoV-2 3CLpro is illustrated. Then, it is shown that the 3C-like proteases of family Coronaviridae (the protease family C30), which are closely related to SARS-CoV-2 3CLpro, have the same homodimeric structure and common activation mechanism via domain III mediated dimerization. The survey extended to order Nidovirales reveals that all 3C-like proteases belonging to Nidovirales have domain III, but with various chain lengths, and 3CLpro of family Mesoniviridae (family C107) has the same homodimeric structure as that of C30, even though they have no sequence similarity. As a reference, monomeric 3C proteases belonging to the more distant family Picornaviridae (family C3) lacking domain III are compared with C30, and it is shown that the 3C proteases are rigid enough to maintain their structures in the active state. Supplementary Information: The online version contains supplementary material available at 10.1007/s12551-022-01020-x.

Multiscale Enhanced Sampling Using Machine Learning.

Multiscale enhanced sampling (MSES) allows for an enhanced sampling of all-atom protein structures by coupling with the accelerated dynamics of the associated coarse-grained (CG) model. In this paper, we propose an MSES extension to replace the CG model with the dynamics on the reduced subspace generated by a machine learning approach, the variational autoencoder (VAE). The molecular dynamic (MD) trajectories of the ribose-binding protein (RBP) in both the closed and open forms were used as the input by extracting the inter-residue distances as the structural features in order to train the VAE model, allowing the encoded latent layer to characterize the difference in the structural dynamics of the closed and open forms. The interpolated data characterizing the RBP structural change in between the closed and open forms were thus efficiently generated in the low-dimensional latent space of the VAE, which was then decoded into the time-series data of the inter-residue distances and was useful for driving the structural sampling at an atomistic resolution via the MSES scheme. The free energy surfaces on the latent space demonstrated the refinement of the generated data that had a single basin into the simulated data containing two closed and open basins, thus illustrating the usefulness of the MD simulation together with the molecular mechanics force field in recovering the correct structural ensemble.

Allosteric Regulation of 3CL Protease of SARS-CoV-2 and SARS-CoV Observed in the Crystal Structure Ensemble.

The 3C-like protease (3CLpro) of SARS-CoV-2 is a potential therapeutic target for COVID-19. Importantly, it has an abundance of structural information solved as a complex with various drug candidate compounds. Collecting these crystal structures (83 Protein Data Bank (PDB) entries) together with those of the highly homologous 3CLpro of SARS-CoV (101 PDB entries), we constructed the crystal structure ensemble of 3CLpro to analyze the dynamic regulation of its catalytic function. The structural dynamics of the 3CLpro dimer observed in the ensemble were characterized by the motions of four separate loops (the C-loop, E-loop, H-loop, and Linker) and the C-terminal domain III on the rigid core of the chymotrypsin fold. Among the four moving loops, the C-loop (also known as the oxyanion binding loop) causes the order (active)-disorder (collapsed) transition, which is regulated cooperatively by five hydrogen bonds made with the surrounding residues. The C-loop, E-loop, and Linker constitute the major ligand binding sites, which consist of a limited variety of binding residues including the substrate binding subsites. Ligand binding causes a ligand size dependent conformational change to the E-loop and Linker, which further stabilize the C-loop via the hydrogen bond between the C-loop and E-loop. The T285A mutation from SARS-CoV 3CLpro to SARS-CoV-2 3CLpro significantly closes the interface of the domain III dimer and allosterically stabilizes the active conformation of the C-loop via hydrogen bonds with Ser1 and Gly2; thus, SARS-CoV-2 3CLpro seems to have increased activity relative to that of SARS-CoV 3CLpro.

Path Ensembles for Pin1-Catalyzed Cis-Trans Isomerization of a Substrate Calculated by Weighted Ensemble Simulations.

Pin1 enzyme protein recognizes specifically phosphorylated serine/threonine (pSer/pThr) and catalyzes the slow interconversion of the peptidyl-prolyl bond between cis and trans forms. Structural dynamics between the cis and trans forms are essential to reveal the underlying molecular mechanism of the catalysis. In this study, we apply the weighted ensemble (WE) simulation method to obtain comprehensive path ensembles for the Pin1-catalyzed isomerization process. Associated rate constants for both cis-to-trans and trans-to-cis isomerization are calculated to be submicroseconds time scales, which are in good agreement with the calculated free energy landscape where the cis form is slightly less favorable. The committor-like analysis indicates the shift of the transition state toward trans form (at the isomerization angle ω ∼ 110°) compared to the intrinsic position for the isolated substrate (ω ∼ 90°). The calculated structural ensemble clarifies a role of both the dual-histidine motif, His59/His157, and the basic residues, Lys63/Arg68/Arg69, to anchor both sides of the peptidyl-prolyl bond, the aromatic ring in Pro, and the phosphate in pSer, respectively. The rotation of the torsion angle is found to be facilitated by relaying the hydrogen-bond partner of the main-chain oxygen in pSer from Cys113 in the cis form to Arg68 in the trans form, through Ser154 at the transition state, which is really the cause of the shift in the transition state. The role of Ser154 as a driving force of the isomerization is confirmed by additional WE and free energy calculations for S154A mutant where the isomerization takes place slightly slower and the free energy barrier increases through the mutation. The present study shows the usefulness of the WE simulation for substantial path samplings between the reactant and product states, unraveling the molecular mechanism of the enzyme catalysis.

Molecular basis of ubiquitin-specific protease 8 autoinhibition by the WW-like domain.

Ubiquitin-specific protease 8 (USP8) is a deubiquitinating enzyme involved in multiple membrane trafficking pathways. The enzyme activity is inhibited by binding to 14-3-3 proteins. Mutations in the 14-3-3-binding motif in USP8 are related to Cushing's disease. However, the molecular basis of USP8 activity regulation remains unclear. This study identified amino acids 645-684 of USP8 as an autoinhibitory region, which might interact with the catalytic USP domain, as per the results of pull-down and single-molecule FRET assays performed in this study. In silico modelling indicated that the region forms a WW-like domain structure, plugs the catalytic cleft, and narrows the entrance to the ubiquitin-binding pocket. Furthermore, 14-3-3 inhibited USP8 activity partly by enhancing the interaction between the WW-like and USP domains. These findings provide the molecular basis of USP8 autoinhibition via the WW-like domain. Moreover, they suggest that the release of autoinhibition may underlie Cushing's disease due to USP8 mutations.

Scalable free energy calculation of proteins via multiscale essential sampling.

A multiscale simulation method, "multiscale essential sampling (MSES)," is proposed for calculating free energy surface of proteins in a sizable dimensional space with good scalability. In MSES, the configurational sampling of a full-dimensional model is enhanced by coupling with the accelerated dynamics of the essential degrees of freedom. Applying the Hamiltonian exchange method to MSES can remove the biasing potential from the coupling term, deriving the free energy surface of the essential degrees of freedom. The form of the coupling term ensures good scalability in the Hamiltonian exchange. As a test application, the free energy surface of the folding process of a miniprotein, chignolin, was calculated in the continuum solvent model. Results agreed with the free energy surface derived from the multicanonical simulation. Significantly improved scalability with the MSES method was clearly shown in the free energy calculation of chignolin in explicit solvent, which was achieved without increasing the number of replicas in the Hamiltonian exchange.