Verifying trivial parallelization of multicanonical molecular dynamics for conformational sampling of a polypeptide in explicit water

Molecular Mechanisms of Functional Modulation of Transcriptional Coactivator PC4 via Phosphorylation on Its Intrinsically Disordered Region

To investigate the effects of phosphorylation on the function of the human positive cofactor 4 (PC4), an enhanced molecular dynamics (MD) simulation was performed. The simulation system consists of the N-terminal intrinsic disordered region (IDR) of PC4 and a complex that comprises the C-terminal acidic activation domain of a herpes simplex virion protein 16 (VP16ad) and a homodimer of the C-terminal structured core domain of PC4 (PC4ctd). An earlier report of an experimental study reported that the PC4-VP16ad interaction is modulated by incremental phosphorylation of the IDR. That report also proposed a dynamic model where phosphorylated serine residues of a segment (SEAC) in the IDR contact positively charged residues (lysin and arginine) of another segment (K-rich) in the IDR. This contact formation induced by the phosphorylation results in variation of PC4-VP16ad interaction. However, this contact formation has not yet been measured directly because it is transiently formed and because the SEAC and K-rich segments are unstructured with high flexibility. We performed two simulations to mimic the incremental phosphorylation: the IDR was not phosphorylated in one simulation and only partially phosphorylated in the other. Our simulation results indicate that the phosphorylation weakens the IDR-VP16ad contact considerably with the induction of a compact structure in the IDR. This structure was stabilized by electrostatic interactions between the phosphorylated serine residues of a segment and lysine or arginine residues of another segment in the IDR, but the conformational fluctuation of this compact structure was considerably large. Consequently, the present study supports the experimentally proposed dynamic model. Results of this study can be important for computational elucidation of the functional modulation of PC4.

Computer simulation of molecular recognition in biomolecular system: from in silico screening to generalized ensembles

Prediction of ligand-receptor complex structure is important in both the basic science and the industry such as drug discovery. We report various computation molecular docking methods: fundamental in silico (virtual) screening, ensemble docking, enhanced sampling (generalized ensemble) methods, and other methods to improve the accuracy of the complex structure. We explain not only the merits of these methods but also their limits of application and discuss some interaction terms which are not considered in the in silico methods. In silico screening and ensemble docking are useful when one focuses on obtaining the native complex structure (the most thermodynamically stable complex). Generalized ensemble method provides a free-energy landscape, which shows the distribution of the most stable complex structure and semi-stable ones in a conformational space. Also, barriers separating those stable structures are identified. A researcher should select one of the methods according to the research aim and depending on complexity of the molecular system to be studied.

Fly casting with ligand sliding and orientational selection supporting complex formation of a GPCR and a middle sized flexible molecule

A GA-guided multidimensional virtual-system coupled molecular dynamics (GA-mD-VcMD) simulation was conducted to elucidate binding mechanisms of a middle-sized flexible molecule, bosentan, to a GPCR protein, human endothelin receptor type B (hETB). GA-mD-VcMD is a generalized ensemble method that produces a free-energy landscape of the ligand-receptor binding by searching large-scale motions accompanied with stable maintenance of the fragile cell-membrane structure. All molecular components (bosentan, hETB, membrane, and solvent) were represented with an all-atom model. Then sampling was conducted from conformations where bosentan was distant from the binding site in the hETB binding pocket. The deepest basin in the resultant free-energy landscape was assigned to native-like complex conformation. The following binding mechanism was inferred. First, bosentan fluctuating randomly in solution is captured using a tip region of the flexible N-terminal tail of hETB via nonspecific attractive interactions (fly casting). Bosentan then slides occasionally from the tip to the root of the N-terminal tail (ligand–sliding). During this sliding, bosentan passes the gate of the binding pocket from outside to inside of the pocket with an accompanying rapid reduction of the molecular orientational variety of bosentan (orientational selection). Last, in the pocket, ligand–receptor attractive native contacts are formed. Eventually, the native-like complex is completed. The bosentan-captured conformations by the tip-region and root-region of the N-terminal tail correspond to two basins in the free-energy landscape. The ligand-sliding corresponds to overcoming of a free-energy barrier between the basins.

Cyclosporin A: Conformational Complexity and Chameleonicity

The chameleonic behavior of cyclosporin A (CsA) was investigated through conformational ensembles employing multicanonical molecular dynamics simulations that could sample the cis and trans isomers of N-methylated amino acids; these assessments were conducted in explicit water, dimethyl sulfoxide, acetonitrile, methanol, chloroform, cyclohexane (CHX), and n-hexane (HEX) using AMBER ff03, AMBER10:EHT, AMBER12:EHT, and AMBER14:EHT force fields. The conformational details were discussed employing the free-energy landscapes (FELs) at T = 300 K; it was observed that the experimentally determined structures of CsA were only a part of the conformational space. Comparing the ROESY measurements in CHX-d12 and HEX-d14, the major conformations in those apolar solvents were essentially the same as that in CDCl₃ except for the observation of some sidechain rotamers. The effects of the metal ions on the conformations, including the cis/trans isomerization, were also investigated. Based on the analysis of FELs, it was concluded that the AMBER ff03 force field best described the experimentally derived conformations, indicating that CsA intrinsically formed membrane-permeable conformations and that the metal ions might be the key to the cis/trans isomerization of N-methylated amino acids before binding a partner protein.

Generalized-ensemble method study: A helix-mimetic compound inhibits protein-protein interaction by long-range and short-range intermolecular interactions

A heterocyclic compound mS-11 is a helix-mimetic designed to inhibit binding of an intrinsic disordered protein neural restrictive silence factor/repressor element 1 silencing factor (NRSF/REST) to a receptor protein mSin3B. We apply a generalized ensemble method, multi-dimensional virtual-system coupled molecular dynamics developed by ourselves recently, to a system consisting of mS-11 and mSin3B, and obtain a thermally equilibrated distribution, which is comprised of the bound and unbound states extensively. The lowest free-energy position of mS-11 coincides with the NRSF/REST position in the experimentally-determined NRSF/REST-mSin3B complex. Importantly, the molecular orientation of mS-11 is ordering in a wide region around mSin3B. The resultant binding scenario is: When mS-11 is distant from the binding site of mSin3B, mS-11 descends the free-energy slope toward the binding site maintaining the molecular orientation to be advantageous for binding. Then, finally a long and flexible hydrophobic sidechain of mS-11 fits into the binding site, which is the lowest-free-energy complex structure inhibiting NRSF/REST binding to mSin3B.

Difference of binding modes among three ligands to a receptor mSin3B corresponding to their inhibitory activities

A preceding experiment suggested that a compound, which inhibits binding of the REST/NRSF segment to the cleft of a receptor protein mSin3B, can be a potential drug candidate to ameliorate many neuropathies. We have recently developed an enhanced conformational sampling method, genetic-algorithm-guided multi-dimensional virtual-system-coupled canonical molecular dynamics, and in the present study, applied it to three systems consisting of mSin3B and one of three compounds, sertraline, YN3, and acitretin. Other preceding experiments showed that only sertraline inhibits the binding of REST/NRSF to mSin3B. The current simulation study produced the spatial distribution of the compounds around mSin3B, and showed that sertraline and YN3 bound to the cleft of mSin3B with a high propensity, although acitretin did not. Further analyses of the simulation data indicated that only the sertraline-mSin3B complex produced a hydrophobic core similar to that observed in the molecular interface of the REST/NRSF-mSin3B complex: An aromatic ring of sertraline sunk deeply in the mSin3B's cleft forming a hydrophobic core contacting to hydrophobic amino-acid residues located at the bottom of the cleft. The present study proposes a step to design a compound that inhibits competitively the binding of a ligand to its receptor.

Elucidating Solution Structures of Cyclic Peptides Using Molecular Dynamics Simulations

Protein-protein interactions are vital to biological processes, but the shape and size of their interfaces make them hard to target using small molecules. Cyclic peptides have shown promise as protein-protein interaction modulators, as they can bind protein surfaces with high affinity and specificity. Dozens of cyclic peptides are already FDA approved, and many more are in various stages of development as immunosuppressants, antibiotics, antivirals, or anticancer drugs. However, most cyclic peptide drugs so far have been natural products or derivatives thereof, with de novo design having proven challenging. A key obstacle is structural characterization: cyclic peptides frequently adopt multiple conformations in solution, which are difficult to resolve using techniques like NMR spectroscopy. The lack of solution structural information prevents a thorough understanding of cyclic peptides' sequence-structure-function relationship. Here we review recent development and application of molecular dynamics simulations with enhanced sampling to studying the solution structures of cyclic peptides. We describe novel computational methods capable of sampling cyclic peptides' conformational space and provide examples of computational studies that relate peptides' sequence and structure to biological activity. We demonstrate that molecular dynamics simulations have grown from an explanatory technique to a full-fledged tool for systematic studies at the forefront of cyclic peptide therapeutic design.

GA-guided mD-VcMD: A genetic-algorithm-guided method for multi-dimensional virtual-system coupled molecular dynamics

We introduced a conformational sampling method in an earlier report: The multi-dimensional virtual-system coupled molecular dynamics (mD-VcMD) enhances conformational sampling of a biomolecular system by computer simulations. Herein, new sampling method, a subzone-based mD-VcMD, is presented as an extension of mD-VcMD. Then, the subzone-based method is extended further using a genetic algorithm (GA) named the GA-guided mD-VcMD. In these methods, iterative simulation runs are performed to increase the sampled region gradually. The new methods have the following benefits: (1) They are free from a production run: i.e., all snapshots from all iterations are useful for analyses. (2) They are free from fine tuning of a weight function (probability distribution function or potential of mean force). (3) A canonical ensemble (i.e., a thermally equilibrated ensemble) is generated from a simple procedure. A thermodynamic weight is assigned to each snapshot. (4) Selective sampling can be performed for particularly addressing a poorly sampled region without breaking the proportion of the canonical ensemble if the poorly sampled conformational region emerges in sampling. By applying the methods to a simple system that involves an energy barrier between potential-energy minima, we demonstrated that the new methods have considerably higher sampling efficiency than the original mD-VcMD does.

Molecular Interaction Mechanism of a 14-3-3 Protein with a Phosphorylated Peptide Elucidated by Enhanced Conformational Sampling

Enhanced conformational sampling, a genetic-algorithm-guided multidimensional virtual-system coupled molecular dynamics, can provide equilibrated conformational distributions of a receptor protein and a flexible ligand at room temperature. The distributions provide not only the most stable but also semistable complex structures and propose a ligand-receptor binding process. This method was applied to a system consisting of a receptor protein, 14-3-3ε, and a flexible peptide, phosphorylated myeloid leukemia factor 1 (pMLF1). The results present comprehensive binding pathways of pMLF1 to 14-3-3ε. We identified four thermodynamically stable clusters of MLF1 on the 14-3-3ε surface and free-energy barriers among some clusters. The most stable cluster includes two high-density spots connected by a narrow corridor. When pMLF1 passes the corridor, a salt-bridge relay (switching) related to the phosphorylated residue of pMLF1 occurs. Conformations in one high-density spot are similar to the experimentally determined complex structure. Three-dimensional distributions of residues in the intermolecular interface rationally explain the binding constant changes resulting from the alanine mutation experiment for the residues. We also performed a simulation of nonphosphorylated peptide and 14-3-3ε, which demonstrated that the complex structure was unstable, suggesting that phosphorylation of the peptide is crucially important for binding to 14-3-3ε.

Free-energy landscape of molecular interactions between endothelin 1 and human endothelin type B receptor: fly-casting mechanism

The free-energy landscape of interaction between a medium-sized peptide, endothelin 1 (ET1), and its receptor, human endothelin type B receptor (hETB), was computed using multidimensional virtual-system coupled molecular dynamics, which controls the system's motions by introducing multiple reaction coordinates. The hETB embedded in lipid bilayer was immersed in explicit solvent. All molecules were expressed as all-atom models. The resultant free-energy landscape had five ranges with decreasing ET1-hETB distance: completely dissociative, outside-gate, gate, binding pocket, and genuine-bound ranges. In the completely dissociative range, no ET1-hETB interaction appeared. In the outside-gate range, an ET1-hETB attractive interaction was the fly-casting mechanism. In the gate range, the ET1 orientational variety decreased rapidly. In the binding pocket range, ET1 was in a narrow pathway with a steep free-energy slope. In the genuine-bound range, ET1 was in a stable free-energy basin. A G-protein-coupled receptor (GPCR) might capture its ligand from a distant place.

Multidimensional virtual-system coupled canonical molecular dynamics to compute free-energy landscapes of peptide multimer assembly

An enhanced-sampling method termed multidimensional virtual-system coupled canonical molecular dynamics (mD-VcMD) method is developed. In many cases, generalized-ensemble methods realizing enhanced sampling, for example, adaptive umbrella sampling, apply an effective potential, which is derived from temporarily assumed canonical distribution as a function of one or more arbitrarily defined reaction coordinates. However, it is not straightforward to estimate the appropriate canonical distribution, especially for cases applying multiple reaction coordinates. The current method, mD-VcMD, does not rely on the form of the canonical distribution. Therefore, it is practically useful to explore a high-dimensional reaction-coordinate space. In this article, formulation of mD-VcMD and its evaluation with the simple molecular models consisting of three or four alanine peptides are presented. We confirmed that mD-VcMD efficiently searched 2D and 3D reaction-coordinate spaces defined as interpeptide distances. Direct comparisons with results of long-term canonical MD simulations revealed that mD-VcMD produces correct canonical ensembles. © 2019 Wiley Periodicals, Inc.

Conformation and Permeability: Cyclic Hexapeptide Diastereomers

Conformational ensembles of eight cyclic hexapeptide diastereomers in explicit cyclohexane, chloroform, and water were analyzed by multicanonical molecular dynamics (McMD) simulations. Free-energy landscapes (FELs) for each compound and solvent were obtained from the molecular shapes and principal component analysis at T = 300 K; detailed analysis of the conformational ensembles and flexibility of the FELs revealed that permeable compounds have different structural profiles even for a single stereoisomeric change. The average solvent-accessible surface area (SASA) in cyclohexane showed excellent correlation with the cell permeability, whereas this correlation was weaker in chloroform. The average SASA in water correlated with the aqueous solubility. The average polar surface area did not correlate with cell permeability in these solvents. A possible strategy for designing permeable cyclic peptides from FELs obtained from McMD simulations is proposed.

Multi-dimensional virtual system introduced to enhance canonical sampling

When an important process of a molecular system occurs via a combination of two or more rare events, which occur almost independently to one another, computational sampling for the important process is difficult. Here, to sample such a process effectively, we developed a new method, named the "multi-dimensional Virtual-system coupled Monte Carlo (multi-dimensional-VcMC)" method, where the system interacts with a virtual system expressed by two or more virtual coordinates. Each virtual coordinate controls sampling along a reaction coordinate. By setting multiple reaction coordinates to be related to the corresponding rare events, sampling of the important process can be enhanced. An advantage of multi-dimensional-VcMC is its simplicity: Namely, the conformation moves widely in the multi-dimensional reaction coordinate space without knowledge of canonical distribution functions of the system. To examine the effectiveness of the algorithm, we introduced a toy model where two molecules (receptor and its ligand) bind and unbind to each other. The receptor has a deep binding pocket, to which the ligand enters for binding. Furthermore, a gate is set at the entrance of the pocket, and the gate is usually closed. Thus, the molecular binding takes place via the two events: ligand approach to the pocket and gate opening. In two-dimensional (2D)-VcMC, the two molecules exhibited repeated binding and unbinding, and an equilibrated distribution was obtained as expected. A conventional canonical simulation, which was 200 times longer than 2D-VcMC, failed in sampling the binding/unbinding effectively. The current method is applicable to various biological systems.

Molecular dynamics coupled with a virtual system for effective conformational sampling

An enhanced conformational sampling method is proposed: virtual-system coupled canonical molecular dynamics (VcMD). Although VcMD enhances sampling along a reaction coordinate, this method is free from estimation of a canonical distribution function along the reaction coordinate. This method introduces a virtual system that does not necessarily obey a physical law. To enhance sampling the virtual system couples with a molecular system to be studied. Resultant snapshots produce a canonical ensemble. This method was applied to a system consisting of two short peptides in an explicit solvent. Conventional molecular dynamics simulation, which is ten times longer than VcMD, was performed along with adaptive umbrella sampling. Free-energy landscapes computed from the three simulations mutually converged well. The VcMD provided quicker association/dissociation motions of peptides than the conventional molecular dynamics did. The VcMD method is applicable to various complicated systems because of its methodological simplicity. © 2018 Wiley Periodicals, Inc.

Dynamics of the Extended String-Like Interaction of TFIIE with the p62 Subunit of TFIIH

General transcription factor II E (TFIIE) contains an acid-rich region (residues 378-393) in its α-subunit, comprising 13 acidic and two hydrophobic (Phe387 and Val390) residues. Upon binding to the p62 subunit of TFIIH, the acidic region adopts an extended string-like structure on the basic groove of the pleckstrin homology domain (PHD) of p62, and inserts Phe387 and Val390 into two shallow pockets in the groove. Here, we have examined the dynamics of this interaction by NMR and molecular dynamics (MD) simulations. Although alanine substitution of Phe387 and/or Val390 greatly reduced binding to PHD, the binding mode of the mutants was similar to that of the wild-type, as judged by the chemical-shift changes of the PHD. NMR relaxation dispersion profiles of the interaction exhibited large amplitudes for residues in the C-terminal half-string in the acidic region (Phe387, Glu388, Val390, Ala391, and Asp392), indicating a two-site binding mode: one corresponding to the final complex structure, and one to an off-pathway minor complex. To probe the off-pathway complex structure, an atomically detailed free-energy landscape of the binding mode was computed by all-atom multicanonical MD. The most thermodynamically stable cluster corresponded to the final complex structure. One of the next stable clusters was the off-pathway structure cluster, showing the reversed orientation of the C-terminal half-string on the PHD groove, as compared with the final structure. MD calculations elucidated that the C-terminal half-acidic-string forms encounter complexes mainly around the positive groove region with nearly two different orientations of the string, parallel and antiparallel to the final structure. Interestingly, the most encountered complexes exhibit a parallel-like orientation, suggesting that the string has a tendency to bind around the groove in the proper orientation with the aid of Phe387 and/or Val390 to proceed smoothly to the final complex structure.

Enhanced conformational sampling to visualize a free-energy landscape of protein complex formation

We introduce various, recently developed, generalized ensemble methods, which are useful to sample various molecular configurations emerging in the process of protein-protein or protein-ligand binding. The methods introduced here are those that have been or will be applied to biomolecular binding, where the biomolecules are treated as flexible molecules expressed by an all-atom model in an explicit solvent. Sampling produces an ensemble of conformations (snapshots) that are thermodynamically probable at room temperature. Then, projection of those conformations to an abstract low-dimensional space generates a free-energy landscape. As an example, we show a landscape of homo-dimer formation of an endothelin-1-like molecule computed using a generalized ensemble method. The lowest free-energy cluster at room temperature coincided precisely with the experimentally determined complex structure. Two minor clusters were also found in the landscape, which were largely different from the native complex form. Although those clusters were isolated at room temperature, with rising temperature a pathway emerged linking the lowest and second-lowest free-energy clusters, and a further temperature increment connected all the clusters. This exemplifies that the generalized ensemble method is a powerful tool for computing the free-energy landscape, by which one can discuss the thermodynamic stability of clusters and the temperature dependence of the cluster networks.

Revisiting antibody modeling assessment for CDR-H3 loop

The antigen-binding site of antibodies, also known as complementarity-determining region (CDR), has hypervariable sequence properties. In particular, the third CDR loop of the heavy chain, CDR-H3, has such variability in its sequence, length, and conformation that ordinary modeling techniques cannot build a high-quality structure. At Stage 2 of the Second Antibody Modeling Assessment (AMA-II) held in 2013, the model structures of the CDR-H3 loops were submitted by the seven modelers and were critically assessed. After our participation in AMA-II, we rebuilt one of the long CDR-H3 loops with 13 residues (A52 antibody) by a more precise method, using enhanced conformational sampling with the explicit water model, as compared to our previous method employed at AMA-II. The current stable models obtained from the free energy landscape at 300 K include structures similar to the X-ray crystal structures. Those models were not built in our previous work at AMA-II. The current free energy landscape suggested that the CDR-H3 loop structures in the crystal are not stable in solution, but they are stabilized by the crystal packing effect.

Enhanced sampling simulations to construct free-energy landscape of protein-partner substrate interaction

Molecular dynamics (MD) simulations using all-atom and explicit solvent models provide valuable information on the detailed behavior of protein-partner substrate binding at the atomic level. As the power of computational resources increase, MD simulations are being used more widely and easily. However, it is still difficult to investigate the thermodynamic properties of protein-partner substrate binding and protein folding with conventional MD simulations. Enhanced sampling methods have been developed to sample conformations that reflect equilibrium conditions in a more efficient manner than conventional MD simulations, thereby allowing the construction of accurate free-energy landscapes. In this review, we discuss these enhanced sampling methods using a series of case-by-case examples. In particular, we review enhanced sampling methods conforming to trivial trajectory parallelization, virtual-system coupled multicanonical MD, and adaptive lambda square dynamics. These methods have been recently developed based on the existing method of multicanonical MD simulation. Their applications are reviewed with an emphasis on describing their practical implementation. In our concluding remarks we explore extensions of the enhanced sampling methods that may allow for even more efficient sampling.

Adaptive lambda square dynamics simulation: an efficient conformational sampling method for biomolecules

A novel, efficient sampling method for biomolecules is proposed. The partial multicanonical molecular dynamics (McMD) was recently developed as a method that improved generalized ensemble (GE) methods to focus sampling only on a part of a system (GEPS); however, it was not tested well. We found that partial McMD did not work well for polylysine decapeptide and gave significantly worse sampling efficiency than a conventional GE. Herein, we elucidate the fundamental reason for this and propose a novel GEPS, adaptive lambda square dynamics (ALSD), which can resolve the problem faced when using partial McMD. We demonstrate that ALSD greatly increases the sampling efficiency over a conventional GE. We believe that ALSD is an effective method and is applicable to the conformational sampling of larger and more complicated biomolecule systems.

Enhanced and effective conformational sampling of protein molecular systems for their free energy landscapes

Protein folding and protein-ligand docking have long persisted as important subjects in biophysics. Using multicanonical molecular dynamics (McMD) simulations with realistic expressions, i.e., all-atom protein models and an explicit solvent, free-energy landscapes have been computed for several systems, such as the folding of peptides/proteins composed of a few amino acids up to nearly 60 amino-acid residues, protein-ligand interactions, and coupled folding and binding of intrinsically disordered proteins. Recent progress in conformational sampling and its applications to biophysical systems are reviewed in this report, including descriptions of several outstanding studies. In addition, an algorithm and detailed procedures used for multicanonical sampling are presented along with the methodology of adaptive umbrella sampling. Both methods control the simulation so that low-probability regions along a reaction coordinate are sampled frequently. The reaction coordinate is the potential energy for multicanonical sampling and is a structural identifier for adaptive umbrella sampling. One might imagine that this probability control invariably enhances conformational transitions among distinct stable states, but this study examines the enhanced conformational sampling of a simple system and shows that reasonably well-controlled sampling slows the transitions. This slowing is induced by a rapid change of entropy along the reaction coordinate. We then provide a recipe to speed up the sampling by loosening the rapid change of entropy. Finally, we report all-atom McMD simulation results of various biophysical systems in an explicit solvent.

Conformational Ensembles of an Intrinsically Disordered Protein pKID with and without a KIX Domain in Explicit Solvent Investigated by All-Atom Multicanonical Molecular Dynamics

The phosphorylated kinase-inducible activation domain (pKID) adopts a helix–loop–helix structure upon binding to its partner KIX, although it is unstructured in the unbound state. The N-terminal and C-terminal regions of pKID, which adopt helices in the complex, are called, respectively, α_A and α_B. We performed all-atom multicanonical molecular dynamics simulations of pKID with and without KIX in explicit solvents to generate conformational ensembles. Although the unbound pKID was disordered overall, α_A and α_B exhibited a nascent helix propensity; the propensity of α_A was stronger than that of α_B, which agrees with experimental results. In the bound state, the free-energy landscape of α_B involved two low free-energy fractions: native-like and non-native fractions. This result suggests that α_B folds according to the induced-fit mechanism. The α_B-helix direction was well aligned as in the NMR complex structure, although the α_A helix exhibited high flexibility. These results also agree quantitatively with experimental observations. We have detected that the α_B helix can bind to another site of KIX, to which another protein MLL also binds with the adopting helix. Consequently, MLL can facilitate pKID binding to the pKID-binding site by blocking the MLL-binding site. This also supports experimentally obtained results.

Theory for trivial trajectory parallelization of multicanonical molecular dynamics and application to a polypeptide in water

Trivial trajectory parallelization of multicanonical molecular dynamics (TTP-McMD) explores the conformational space of a biological system with multiple short runs of McMD starting from various initial structures. This method simply connects (i.e., trivially parallelizes) the short trajectories and generates a long trajectory. First, we theoretically prove that the simple trajectory connection satisfies a detailed balance automatically. Thus, the resultant long trajectory is regarded as a single multicanonical trajectory. Second, we applied TTP-McMD to an alanine decapeptide with an all-atom model in explicit water to compute a free-energy landscape. The theory imposes two requirements on the multiple trajectories. We have demonstrated that TTP-McMD naturally satisfies the requirements. The TTP-McMD produces the free-energy landscape considerably faster than a single-run McMD does. We quantitatively showed that the accuracy of the computed landscape increases with increasing the number of multiple runs. Generally, the free-energy landscape of a large biological system is unknown a priori. The current method is suitable for conformational sampling of such a large system to reduce the waiting time to obtain a canonical ensemble statistically reliable.