Folding simulation of small proteins by dissipative particle dynamics (DPD) with non-empirical interaction parameters based on fragment molecular orbital calculations

Recently, we have developed a new simulation scheme with dissipative particle dynamics (DPD) based on non-empirical interaction parameters derived from a series of fragment molecular orbital (FMO) calculations. This approach (FMO–DPD) was applied to folding simulations of Chignolin and Superchignolin mini-proteins. Their characteristic hairpin structure was obtained from an elongated form within short computation time. Essential residue-residue interactions such as hydrogen bonding and CH/π were observed at the final form. FMO–DPD should have a potential applicability to nano-bio systems involving proteins.

P roteins are fundamental substances widely involving in biological life phenomena, and their functionalities depend on the folded 3-dimensional structures. Nanobio technology devices with sensor or receptor proteins have attracted practical interests in the field of applied physics as well. In molecular simulations, the reproduction of folded protein structures has long been of the most demanding problem, because of large-scale configurational changes should take place. Amounts of studies with all-atom molecular dynamics (MD) have thus been made. In particular, Shaw's group 1) provided a revolutionary approach with customized computer chip (ANTON) to accelerate MD simulation and succeeded in folding simulations as long as ms or longer. The folding simulations with atomistic MD are still costly, unless using supercomputers or special accelerators. Therefore, the so-called coarse-grained (CG) approaches have attracted interests as alternative routes. For example, Trp-Cage 2) (consisting of 20 amino acid residues) was successfully folded by CG-MD, 3) and the method of dissipative particle dynamics (DPD) 4) was applied to handle the folding of neurotoxic α-synuclein 5) (140 residues). DPD could have a high efficiency in surveying wide configuration spaces of protein structures (not limited to the folding). As a brand new method, the folding of model proteins was solved through quantum annealing. 6,7) Although the CG simulations for proteins were certainly cost-effective, the sets of parameters describing various interactions (among residue-residue and residue-water) were of empirical type. In other words, the applicability and usability could be limited potentially.
Recently, we have developed a non-empirical simulation scheme of DPD in which the set of effective interaction parameters a ij associated with χ-parameters are evaluated by the fragment molecular orbital (FMO) calculations, 8) abbreviated to FMO-DPD. The corresponding workflow to obtain χ-parameters by FMO was available as a generic system named FCEWS (FMO-based Chi-parameter Evaluation Workflow System). 9,10) Additionally, we have made an original DPD program, CAMUS 11,12) by which specific interactions such as 1-3 hydrogen bonds 5) are easily incorporated. In this work, the FMO-DPD simulation was applied to the folding of Chignolin 13) composed of 10 amino acid residues. Very recently, detailed dynamics of Chignolin has been reported. 14) The speed-up of such extensive MD simulations still remains an important issue, however. Our FMO-DPD may be an efficient CG approach to the folding problem.
We performed the FMO-DPD simulation of Chignolin as follows. First, the amino acid residues of Chignolin (Gly-Tyr-Asp-Pro-Glu-Thr-Gly-Thr-Trp-Gly) were segmented to the respective particles in DPD (Fig. 1). According to the recipe of Ref. 5, the main-chain backbone was described in Gly unit, and the side chain moieties as well as the N-and Cterminals were individually treated with the respective chemical structures. The non-bonding interaction parameters among the residue particles and the environmental water particle were determined through the FMO-based χ-parameter evaluations. 9,10,15) The segment structures were preoptimized with the GAUSSIAN 09 program 16) at the B97D/ 6-31G(d′) level. The combination of FCEWS 9,10) and ABINIT-MP (our FMO program) 17) was employed in a series of evaluations of χ-parameters, where the interaction energies among segmented particles were computed at the second-order perturbation [or MP2/6-31G(d′)] level. The hydration effect was incorporated by the Poisson-Boltzmann model (FMO-PB). 18,19) In the DPD simulations with the CAMUS code, 11) the force of ith particle f i contains four parts as The first three forces of the right-hand-side are considered within a certain cutoff radius r . c The conservative force F ij C is a soft repulsion action as follows Content from this work may be used under the terms of the Creative Commons Attribution 4.0 license. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
where a ij is the maximum repulsion force between particles i and j, and related definitions are =r r r, The repulsion parameters a ij between particles of different types are related to the Flory-Huggins c-parameter as 4) where a a ii jj ( ) describes the self-interaction of particle i ( j). In Eq. (1), the dissipative force F ij D and the random force F ij R represent hydrodynamic drags and thermal noises of the Gaussian statistics, respectively. The fourth term in the righthand side of Eq. (1) provides an additional spring force F ij S for directly bonded particles in polymers. For a certain connected particle pair i and j, the corresponding Harmonic force F ij H and Morse force F ij M are given, respectively, as In our simulation, the a ij parameters obtained from FMO calculations were used as essential interactions between nonbonded particles (refer to Table I). Note here that Trp (W) side chain was represented with a couple of particles because of its relative largeness (see Fig. 1). The 1-2 bonds and 1-3 bonds were of Harmonic type with C = 160 kT, r e = 0.6 DPD-length unit (R c ) and C = 80 kT, r e = 1.2 R c , respectively. 5) To maintain the relative positions between main-and side-chains, a Morse potential was introduced for the 1-3 bonds (between main-chain Gly (G) particle and neighbored side chain particles) with setting of K M = 12, α = 8, and r e = 0.85 R c (our original setting). Similarly, another Morse potential (K M = 12, α = 8, and r e = 1.34 R c ) was used for the 1-4 bonds. The DPD time step was set as Δt = 0.005 DPD-time unit (at least 0.7 ps in real unit), and the number of particles in a cell [6.9 R c (4.9 nm) size box cell] was 1000. The initial structure of mini-protein was an elongated line structure.
The numerical results in the a ij form are listed for the combinations of particles in Table I, where two values for Gly (N-terminal) is set to zero because the calculated value was negative (or too attractive). Taking the C-terminal end of Gly (G) as an example, it is shown that the interaction value with itself (a ii ) shows a large value of repulsion with 145 kT and that the similar repulsions in a ij are found for the negatively charged Asp (E) and Glu (D). Interestingly, the values involving water took various values, reflecting the respective chemical characteristics of the interacting particles. 15) The bond potentials of 1-2 (direct bond), 1-3 (angle), and 1-4 (dihedral) were appropriately set to maintain the rigidity of protein skeleton in CAMUS. 11,12) Figure 2(a) illustrates the hairpin-shape structure of folded Chignolin, 13) where the internal attractive interactions among residues such as hydrogen bond between Asp3 (E) and Gly7 (G) and CH/π interaction between Tyr2 (Y) and Trp9 (W) are crucial. Several DPD simulations with the FMO-derived effective parameters provided the folded structures, as shown in Figs. 2(b) and 2(c), starting from the elongated ones. It is notable that the computation time for this folding simulation was only 8 s (running 10000 steps) on 8 cores (8 OpenMP threads of parallelization 11,12) ) of an in-house Xeon (E5-2690) server. The FMO-DPD simulation was applied also to folding of Superchignolin 20) (Tyr-Tyr-Asp-Pro-Glu-Thr-Gly-Thr-Trp-Tyr), and the hairpin like structure was reproduced successfully, as illustrated in Fig. 3. Figure 4 plots the distances between N-and C-terminal particles against the time step of DPD for both Chignolin and Superchignolin, indicating the rapid decrease of distance from the elongated structure.
Comparison between Chignolin and Superchignolin suggests that the step-dependent fluctuation of Superchignolin is smaller than that of Chignolin, and this situation is consistent with the fact that the internal interactions in the former is stronger (by the intended design 20) ) than that in the latter.     In this paragraph, verification of geometrical features obtained by our DPD simulation is described for Chignolin. In comparison with the NMR structure, 13) Fig. 5 plots the root mean square deviation (RMSD) of Cα atoms for 10 trajectories (up to 100 000 steps), where the average value (bold blue color plot) reaches the range of 3-4 Å at about 5000 steps and the minimal value is evaluated as 1.6 Å (in black color plot). These values may be acceptable, because both interaction parameters and distances between particles were not set at atomistic level in DPD (unlike usual MD). Related verification with RMSD plots was done as well, and the results are given in the supporting information (Figs. S1 and S2), available online at stacks.iop.org/APEX/13/017002/ mmedia. Another comparison is made for time evolutions of mutual positions (see the corresponding Fig. S3) of Asp3 -Gly7 pair (hydrogen bonding) and Tyr2 -Trp9 pair (weaker CH/π interaction). Rather stable behavior is found for the former, but several different configurations (including misfolded type as shown in Fig. S4) are observed for the latter. These characteristics are consistent with atomistic MD results of Ref. 21.
In this brief letter, we have reported an application of the FMO-DPD simulation to the folding of Chignolin and Superchignolin, 1,3) where the non-bonding particle-particle interaction parameters were determined in non-empirical way. 9,10) The folding was accomplished in short DPD runs, indicating an efficiency in handling proteins. FMO-DPD was already successfully applied by us to lipid membranes and vesicles, [22][23][24] and the technique of reverse mapping 25) (by which nano-scale atomic structures are to be reproduced from meso-scale CG structures) has been recently settled for lipids. Preliminary works for lipid membranes containing miniproteins have been in progress. Although nonbonding interaction parameters have been determined with the FMO calculations in the present study, direct bonding parameters as well as structure-constraint parameters have still been of empirical type. 5) Good balance among all these parameters is important in reproducing proper structures (refer to previous Figs. S1 and S2 and also Fig. S5 of α-synuclein). Nonempirical determination of bonding and constraint parameters should be our future issue. We would continue a series of developments associated with FMO-DPD simulations, considering that this protocol should have a potential to efficiently model the crucial parts (e.g. interfaces) of nanobio device systems consisting of lipids, proteins and substrates.