Simulations of Functional Motions of Super Large Biomolecules with a Mixed-Resolution Model

Abstract

Many large protein machines function through an interplay between large-scale movements and intricate conformational changes. Understanding functional motions of these proteins through simulations becomes challenging for both all-atom and coarse-grained (CG) modeling techniques because neither approach alone can readily capture the full details of these motions. In this study, we develop a multiscale model by employing the popular MARTINI CG model to represent a heterogeneous environment and structurally stable proteins and using the united-atom (UA) model PACE to describe proteins undergoing subtle conformational changes. PACE was previously developed to be compatible with the MARTINI solvent and membrane. Here, we couple the protein descriptions of the two models by directly mixing UA and CG interaction parameters to greatly simplify parameter determination. Through extensive validations with diverse protein systems in solution or membrane, we demonstrate that only additional parameter rescaling is needed to enable the resulting model to recover the stability of native structures of proteins under mixed representation. Moreover, we identify the optimal scaling factors that can be applied to various protein systems, rendering the model potentially transferable. To further demonstrate its applicability for realistic systems, we apply the model to a mechanosensitive ion channel Piezo1 that has peripheral arms for sensing membrane tension and a central pore for ion conductance. The model can reproduce the coupling between Piezo1's large-scale arm movement and subtle pore opening in response to membrane stress while consuming much less computational costs than all-atom models. Therefore, our model shows promise for studying functional motions of large protein machines.

Figure 1

PACEm coupling scheme. (A) Schematic representation of the UA and CG partitioning. (B) PACEm representation of four soluble proteins (magenta: UA region; cyan: CG region).

Figure 2

Representation of the Piezo1 PACEm model. The pore region and the three beams (in magenta) were depicted by the PACE UA model. The three arms are depicted by the ELNEDYN CG model (in cyan).

Figure 3

Variation in the structural stability of monomeric proteins in solution with respect to γ. (A) The variations in each type of mean RMSDs with different choices of scaling factor γ. (B) The numbers of constructs that exhibited their least RMSD in simulations at specific γ values. (C) Average counts of backbone contacts between UA and CG subunits over all the constructs examined. A backbone contact occurs for a short distance (less than 6 Å) between a Cα site of a UA subunit and a BB site of a CG subunit. The contacts are normalized with respect to those obtained at γ = 0.7. The error bars indicate the standard errors of the mean.

Figure 4

PACEm simulations of four soluble proteins with scaling factor γ = 0.7. (A) Protein structure at 50 ns of PACEm simulation with γ = 0.7. The UA and CG regions are colored in magenta and cyan, respectively. The experimental structure is represented in a green cartoon format. (B) RMSF of the backbone (BB/Cα) for four soluble proteins with respect to its experimental structure, plotted against the residue index. Black and green curves denote the RMSF results using PACEm and the CHARMM36m AA force field, respectively. The magenta and cyan areas denote the regions modeled at UA and CG resolutions, respectively. The error bars denote the standard errors of the mean estimated from three independent simulations for each protein.

Figure 5

Variation in the structural stability of dimeric protein complexes in solution with respect to γ. (A) The variations in each type of mean RMSDs with different choices of scaling factor γ. (B) The numbers of constructs that exhibited their least RMSD in simulations at specific γ values. (C) Average counts of backbone contacts between UA and CG subunits over all the constructs examined. A backbone contact occurs for a short distance (less than 6 Å) between a Cα site of a UA subunit and a BB site of a CG subunit. The contacts are normalized with respect to those obtained at γ = 0.7. The error bars indicate the standard errors of the mean.

Figure 6

Tuning η and ζ to improve PACEm. (A) The mean RMSD (Å) obtained from simulations of proteins at the UA level in the CG membrane with different (η,ζ) values being used. The mean was averaged over the RMSD values for ARI, DAP12, and KcsA, each of which was averaged over the last halves of three independent simulations. The standard errors of the mean were also provided. (B) Backbone HB counts for ARI (left), DAP12 (middle), and KcsA (right) at different ζ’s. η was set to 0.95. Red dashed lines denote the reference values derived from the CHARMM36m AA simulations.

Figure 7

Modeling of ARI with a mixed-resolution representation in the CG membrane. (A) Snapshot of one of the ARI constructs. The UA part of ARI is shown as purple ribbons, and the CG part is shown as cyan sticks. Lipid molecules are shown as yellow sticks, and CG water solvents are shown as green spheres. (B) The variations in each type of mean RMSDs with different choices of scaling factor γ. The two panels denote the results of the two constructs of ARI. The error bars denote the standard errors of the mean over three independent simulations.

Figure 8

Free energies of association (kcal/mol) for three dimer complexes modeled with PACEm when γ = 0.7. (A) Mixed representation of three protein complexes. The CG monomer is colored cyan, and the UA monomer is colored magenta. (B) The PMFs of complex association. The green dashed line indicates association free energy obtained from experiments. In the results for H-Ras/Raf, the black curve denotes the PMF obtained when the presentation of the complex is construct 1: H-Ras was modeled at the CG resolution, and Raf was modeled at the UA resolution. The red curve denotes the PMF obtained after the representation resolutions of the two domains were swapped (construct 2). For homodimers like insulin dimer and EphA1, only a single PMF can be obtained, which is shown as a black curve.

Figure 9

Piezo1 MD simulations with optimized PACEm. (A) Backbone RMSD of the pore and arms regions across three independent native pressure (+1 bar) simulations and three independent membrane tension (−10 bar) simulations. The experimental structures (PDB ID: 6B3R) were used for the RMSD calculations. (B) The projection area of three independent native pressure (+1 bar) simulations and three independent membrane tension (−10 bar) simulations. The area is calculated using the formula π × r², where r is derived from . Here, a, b, and c represent the distances between the backbone particle from residue ILE859, located within the outermost helix embedded in the bilayer from three Piezo1 arms. The ball marker indicates the position of ILE859. (C) The flattening angle of the three Piezo1 arms across three independent native pressure (+1 bar) simulations and three independent membrane tension (−10 bar) simulations. The angle is defined by the angle between the arm axis (determined by the COM of the outermost helix, residue 850–860, of the arm and the pore region) and the internal axis (determined by the COM of the cap and CTD region), as illustrated on the left. The values of arm 1, arm 2, and arm 3 from the final 100 ns are plotted here. The average value is indicated by the “×” symbol. (D) Snapshots at the 0 and 200 ns of PACEm simulations (replica 1) of Piezo1 under native pressure conditions (top) and membrane tension conditions (bottom), respectively.

Figure 10

Conformational changes of the Piezo1 pore region. (A) Detailed depiction of the pore region and two types of measures of pore sizes based on the positions of pore-lining helices: either the distances (d₁³⁷,d₂³⁷,d₃³⁷) of the three outer helices TM37 to the pore center or the distances (d₁³⁸,d₂³⁸,d₃³⁸) of the three inner helices TM38 to the pore center. Here, the pore center is defined as the centroid of all TM37s and TM38s, and the centroids of individual TM37s and TM38s are used to calculate their distances to the pore enter. (B) Plots of d₃₇ = (d₁³⁷ + d₂³⁷ + d₃³⁷)/3 against d₃₈ = (d₁³⁸ + d₂³⁸ + d₃³⁸)/3 for all the simulations. Each dot corresponds to one snapshot sampled from the simulations, and the dots from different simulations are colored differently. (C) Definition of the pore size based on pore-lining residues V2476: the average distance (d_V = (d₁^V + d₂^V + d₃^V)/3) of the Cβ atoms of the three valine residues to the centroid of these Cβ atoms. The panel on the right denotes the plots of (d₁^V,d₂^V,d₃^V) for all the simulations. The black dots denote the distance values calculated based on the experimental structure of the closed state. (D) The water molecules (shown as red spheres) surrounding the pore region under the native pressure (left) and membrane tension (right). The displayed conformations derive from the first independent replica of the two systems at the 200 ns. (E) The number density of water molecules along the Y-Z panel. The density was accumulated using the last 10 frames of three independent trajectories at a 10 ns interval.