Structure-Preserving Coarse-Grained Simulation of Proteins in Explicit Solvent

Simulation of Protein Structure using a Coarse-Grained Potential incorporating the Backbone Dihedral Interactions | bioRxiv /* */ /* */ <!-- <!-- /*! * yepnope1.5.4 * (c) WTFPL, GPLv2 */ (function(a,b,c){function d(a){return"[object Function]"==o.call(a)}function e(a){return"string"==typeof a}function f(){}function g(a){return!a||"loaded"==a||"complete"==a||"uninitialized"==a}function h(){var a=p.shift();q=1,a?a.t?m(function(){("c"==a.t?B.injectCss:B.injectJs)(a.s,0,a.a,a.x,a.e,1)},0):(a(),h()):q=0}function i(a,c,d,e,f,i,j){function k(b){if(!o&&g(l.readyState)&&(u.r=o=1,!q&&h(),l.onload=l.onreadystatechange=null,b)){"img"!=a&&m(function(){t.removeChild(l)},50);for(var d in y[c])y[c].hasOwnProperty(d)&&y[c][d].onload()}}var j=j||B.errorTimeout,l=b.createElement(a),o=0,r=0,u={t:d,s:c,e:f,a:i,x:j};1===y[c]&&(r=1,y[c]=[]),"object"==a?l.data=c:(l.src=c,l.type=a),l.width=l.height="0",l.onerror=l.onload=l.onreadystatechange=function(){k.call(this,r)},p.splice(e,0,u),"img"!=a&&(r||2===y[c]?(t.insertBefore(l,s?null:n),m(k,j)):y[c].push(l))}function j(a,b,c,d,f){return q=0,b=b||"j",e(a)?i("c"==b?v:u,a,b,this.i++,c,d,f):(p.splice(this.i++,0,a),1==p.length&&h()),this}function k(){var a=B;return a.loader={load:j,i:0},a}var l=b.documentElement,m=a.setTimeout,n=b.getElementsByTagName("script")[0],o={}.toString,p=[],q=0,r="MozAppearance"in l.style,s=r&&!!b.createRange().compareNode,t=s?l:n.parentNode,l=a.opera&&"[object Opera]"==o.call(a.opera),l=!!b.attachEvent&&!l,u=r?"object":l?"script":"img",v=l?"script":u,w=Array.isArray||function(a){return"[object Array]"==o.call(a)},x=[],y={},z={timeout:function(a,b){return b.length&&(a.timeout=b[0]),a}},A,B;B=function(a){function b(a){var a=a.split("!"),b=x.length,c=a.pop(),d=a.length,c={url:c,origUrl:c,prefixes:a},e,f,g;for(f=0;f<d;f++)g=a[f].split("="),(e=z[g.shift()])&&(c=e(c,g));for(f=0;f<b;f++)c=x[f](c);return c}function g(a,e,f,g,h){var i=b(a),j=i.autoCallback;i.url.split(".").pop().split("?").shift(),i.bypass||(e&&(e=d(e)?e:e[a]||e[g]||e[a.split("/").pop().split("?")[0]]),i.instead?i.instead(a,e,f,g,h):(y[i.url]?i.noexec=!0:y[i.url]=1,f.load(i.url,i.forceCSS||!i.forceJS&&"css"==i.url.split(".").pop().split("?").shift()?"c":c,i.noexec,i.attrs,i.timeout),(d(e)||d(j))&&f.load(function(){k(),e&&e(i.origUrl,h,g),j&&j(i.origUrl,h,g),y[i.url]=2})))}function h(a,b){function c(a,c){if(a){if(e(a))c||(j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}),g(a,j,b,0,h);else if(Object(a)===a)for(n in m=function(){var b=0,c;for(c in a)a.hasOwnProperty(c)&&b++;return b}(),a)a.hasOwnProperty(n)&&(!c&&!--m&&(d(j)?j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}:j[n]=function(a){return function(){var b=[].slice.call(arguments);a&&a.apply(this,b),l()}}(k[n])),g(a[n],j,b,n,h))}else!c&&l()}var h=!!a.test,i=a.load||a.both,j=a.callback||f,k=j,l=a.complete||f,m,n;c(h?a.yep:a.nope,!!i),i&&c(i)}var i,j,l=this.yepnope.loader;if(e(a))g(a,0,l,0);else if(w(a))for(i=0;i (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0];var j=d.createElement(s);var dl=l!='dataLayer'?'&l='+l:'';j.src='//www.googletagmanager.com/gtm.js?id='+i+dl;j.type='text/javascript';j.async=true;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-M677548'); Skip to main content Home About Submit ALERTS / RSS Search for this keyword Advanced Search New Results Simulation of Protein Structure using a Coarse-Grained Potential incorporating the Backbone Dihedral Interactions Kanika Kole , View ORCID Profile Abhik Ghosh Moulick , Jaydeb Chakrabarti doi: https://doi.org/10.1101/2025.08.20.671185 Kanika Kole † Department of Physics of Complex Systems, S. N. Bose National Centre for Basic Sciences , Block-JD, Sector-III, Salt Lake, Kolkata-700106, India ‡ Department of Chemistry, Indian Institute of Technology Kharagpur , Kharagpur 721302, India Find this author on Google Scholar Find this author on PubMed Search for this author on this site Abhik Ghosh Moulick † Department of Physics of Complex Systems, S. N. Bose National Centre for Basic Sciences , Block-JD, Sector-III, Salt Lake, Kolkata-700106, India ¶ Institute of Nanotechnology, Karlsruhe Institute of Technology (KIT) , Kaiserstraße 12, 76131 Karlsruhe, Germany Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Abhik Ghosh Moulick Jaydeb Chakrabarti § Department of Chemical and Biological Sciences, S. N. Bose National Centre for Basic Sciences , Block-JD, Sector-III, Salt Lake, Kolkata-700106, India Find this author on Google Scholar Find this author on PubMed Search for this author on this site For correspondence: jaydebchakrabarti{at}gmail.com Abstract Full Text Info/History Metrics Supplementary material Preview PDF Abstract Many biologically relevant processes occur on time and length scales which are far beyond the reach of atomistic simulations. These processes include large protein dynamics and the self-assembly of biological materials. Coarse-grained molecular modeling allows computer simulations on length and time scales 2–3 orders of magnitude larger than atomistic simulations, bridging the gap between the atomistic and meso-scopic scales. However, the structural information involving the dihedral angles is lost in coarse-graining. We develop a simple coarse-grained protein model with structural information in explicit solvent. We represent the center of mass of each residue as a polymer bead and water oxygen as a solvent bead. Each polymer bead has five degrees of freedom: position of the center and two additional variables for the backbone dihedral angles. All interaction parameters for bonded, non-bonded, dihedral coupling and bead-solvent interactions are derived from the equilibrated all-atom molecular dynamics simulation trajectory. We find that our coarse-grained approach reproduces residue-level structural information that closely matches the crystal structures and all-atom simulation results. 1. Introduction Coarse-grained (CG) models describe complex systems at large length and time scales 1 – 7 far beyond the reach of full microscopic all-atom (AA) models. 8 This is achieved by reducing the number of degrees of freedom where a number of atoms are grouped together into a unit. 9 , 10 However, CG models have inherent limitations. Since a number of atoms are grouped together, the explicit information on planes constituted by different atoms are often not retained. Consequently, the angle between atomic planes, also called the dihedral angles, 11 which describe the spatial conformations of complex molecules cannot be adequately described within the CG models. Hence, CG models are largely inadequate to describe phenomena involving bio-macromolecular conformations. The backbone dihedral angles describe protein conformations. There are plenty of cases where the protein conformation lies at the center stage. Protein–ligand interactions, for instance, depend critically on protein conformational states. Proteins which lack a stable, well-defined three-dimensional structure, like the intrinsically disordered proteins (IDP) 12 undergo fluctuations among conformations even in physiological conditions. 13 They often have a tendency to aggregate in the solution phase, 14 leading to pathological conditions. 14 – 16 Sometimes proteins, like milk protein alpha-lactalbumin in denaturing conditions, form a molten globule state in which local stable structure gets disrupted while retaining overall secondary structure. Molten globule proteins are capable of ligand binding 14 , 15 and form aggregation in solution phases. The structural aspects for such proteins are extremely important to understand their function and solution phase properties. On the other hand, due to involvement of several protein molecules in aggregate formation, a full atomistic simulations are almost next to impossible, pointing out the necessity of CG models which can reliably capture information on protein conformation. A wide variety of CG models has been introduced to capture specific aspects for biomolecular systems, like proteins, 17 nucleic acids, 17 ? , 18 lipid membranes, 19 , 20 carbohydrates, 21 water, 22 , 23 and so on. In highly simplified lattice protein-like hydrophobic-hydrophilic (HP) models, 24 – 27 each amino acid is either represented as a solvophobic or polar entity, but ignores protein backbone information. More detailed CG protein models adopt different levels of simplified polypeptide representation. 28 The main chain is represented by either all heavy atoms or one to two united atoms per residue, while the side chain is typically replaced by one or two united atoms. 28 – 30 Realistic CG models typically derive their parameters from detailed theoretical models, such as atomistic or quantum mechanical simulations, 8 following a bottom-up structure based approach as in the SIRAH CG force field. 31 , 32 There are CG models, like the Martini force field 13 , 33 , 34 which adopt a top-down approach, where the model is built in such a way that it can reproduce a set of experimental macroscopic properties. 35 Other well-known CG models for protein simulations, such as SPICA 36 and UNRES, 37 represent alternative strategies those differ in resolution and solvent treatment. SPICA employs an explicit-solvent representation with a small number of beads per residue and combines bottom-up, top-down, and structure-based parameterization, with secondary structure typically stabilized through elastic network restraints. In contrast, UNRES is an implicit-solvent model with a minimal two-bead protein representation that is designed to be foldable and has been applied to a wide range of bio-molecular systems. However, each of these approaches has inherent limitations to obtain reliable dihedral angle as well as secondary-structure information directly. For force fields such as MARTINI or SIRAH, this limitation is commonly addressed through backmapping, 38 a procedure in which atomistic (AA) structures are reconstructed from CG snapshots or trajectories. However, this is quite challenging because back-mapping is not unique and often depends on methodological choices and structural restraints, which can introduce biases and additional sources of uncertainty. There is a clear need for a CG force fields which retain sufficient intrinsic information such as backbone dihedral angle without relying on back-mapping. Goddard III and coworkers introduce a CG Dihedral Probability Grid Monte Carlo (DPG-MC) approach to generate polypeptides 39 and protein structures. 40 C α atoms are used to represent the protein backbone and pairwise model side-chain–side-chain interactions are employed. An internal-coordinate search algorithm is employed that is guided by residue-specific dihedral angle probability distributions derived from known protein structures in the Protein Data Bank (PDB). The initial Cartesian coordinates of the polypeptides and proteins are constructed using BIOGRAF, after which the secondary structures are randomly assigned to each residue. The net charge of the peptide or protein is determined using a Charge-Equilibration scheme. The energy of the resulting structure is then minimized to convergence using the steepest-descent method, followed by conjugate-gradient minimization. The structures are updated using the DPG-MC simulations. At the end of each DPG-MC run, the lowest-energy conformation is further refined by energy minimization using steepest-descent followed by conjugate-gradient methods until convergence. This model thus integrate protein database knowledge with conformational search methodologies. However, the available PDB data are limited, particularly for intrinsically disordered proteins (IDPs). Moreover, since solvent atoms are not treated explicitly, the interactions between the amino-acid residues and solvent cannot be explicitly treated. With this backdrop, we build a simple CG polymer model for protein to capture the structural utilizing statistical mechanics based physical interactions in fully microscopic all-atom (AA) simulations. For a given coarse-grained variable Γ, we compute the distribution of the variable H (Γ) over equilibrated conformations of the protein in the AA simulations. The equilibrium value and energy associated with Γ are calculated using the Boltzmann formula F (Γ) = − k B T lnH (Γ). Fig. 1 shows a schematic of our model, illustrating polymer beads, connected by springs with finite stretching and bending elastic constants and water oxygen by isolated spheres. The features of the model are as follows: (1) Each polymer bead, connected via a chain of finite elastic constants, represents centre of mass of a protein residue as in the primary sequence of the protein. We assign five variables to each bead: three for the coordinate of the center of mass (COM) and two degrees of freedom corresponding to the main chain backbone dihedral angles ( ϕ and ψ ) to account for the backbone structure. (2) Each solvent water molecule is represented explicitly by the oxygen atom. (3) The bead-solvent interactions are calculated from the distribution of water molecules around the COM of a residue over equilibrium AA trajectory. The beads are classified based on solvent interactions: The solvophobic beads repel water oxygen which correspond to the hydrophobic residues, while the sovophilic beads attract water oxygen, corresponding to the hydrophilic residues in a given solvent condition. (4) The dihedral interactions are computed from the joint probability distributions of different dihedral angles over AA trajectory and grouped together depending on solvophobic and solvophilic property of the residues. All other relevant CG model parameters, like the bead size, the backbone elastic constants, bead-bead non-bonded interactions and solvent-solvent interactions are computed using distributions of the relevant variables over equilibrated AA simulation trajectories. Thus, unlike DPG-MC model, our model is not residue specific but depends on their interactions with solvent and takes care of physically relevant interactions. Our main goal is to have a minimal CG interaction model to reproduce the secondary structural elements of the residues reliably. We then perform MC simulations of the model polymer based on Metropolis sampling to generate equilibrium structure of the chain 41 , 42 in which Cartesian coordinates and backbone dihedral angles are updated at each MC step to generate conformations as per the energy costs in the CG model interactions. Download figure Open in new tab Figure 1: A schematic illustrating the main-chain dihedral angles ϕ and ψ associated with each polymer bead. We test our method on GB3, a small protein with 56 amino acid residues, widely used as a model structure 43 in CG simulations due to its small size, well-defined structure and well-documented folding mechanism. We calculate CG model parameters from the AA simulations data of GB3. We then perform CG MC simulations using the CG parameters determined from the AA trajectory of GB3 and find that the structural elements of the proteins are well captured when compared to their crystal structure and AA data. We further apply our method to other well-structured proteins where the CG parameters are taken from GB3 AA trajectory: (1) Homeodomain: a small, 58-residue protein that binds to specific DNA sequences and functions as a transcription factor, playing a crucial role in gene regulation, (2) Ubiquitin, a 76 residue globular protein involved in the targeted degradation of cyclins and other regulatory proteins and (3) Di-ubiquitin (linked via K48), a two-domain protein consisting of 148 residues that functions as a signal for proteasomal degradation. We find that the CG model parameter for GB3 describes well the structures of these folded proteins. Next, we apply our methodology to intrinsically disordered proteins (IDPs). We study a couple of IDPs, α -synuclein ( α S) and λ N. α S is a small abundant neuronal protein in the brain, best known for its role in synaptic function and pathological aggregation in Parkinson’s disease. λ N is an anti-termination protein that recognizes the nut site on the phage RNA and enables transcription to continue past termination signals. AA simulations of α S yield the elastic constants which are smaller compared to well-folded proteins reflecting their high flexibilities. The CG parameters obtained for α S show primarily unfolded structures in both α S and λ N in agreement to the AA data. 2. Methods 2.1. System preparation We study the following systems in both AA and CG simulations: (1) GB3 (PDB ID: 2OED), 44 a small globular protein consisting of 56 residues, (2) homeodomain protein, another small globular protein containing 58 residues, the initial structure of homeodomain protein ( α 2D) is taken from the PDB structure 1K61 45 without DNA and another home-odomain protein ( α 2B), (3) ubiquitin (PDB ID: 1UBQ), 46 a small structural protein containing 76 residues, (4) di-ubiquitin (PDB ID: 7S6O, linked via K48), 47 a multidomain protein consisting of 148 residues, (5) α -synuclein (PDB ID: 1XQ8), 48 an IDP, where only considered the N-terminal (residues 1–60) and NAC regions (residues 61–95), the C-terminal segment (residues 96–140), which is already disorder in the crystal structure is excluded and (6) λ N protein (PDB ID: 1QFQ), 49 an IDP protein containing 35 residues, the initial model structure of this disordered Bacteriophage λ N protein is taken from the RNA-unbound form of the PDB structure 1QFQ. 49 λ N adopts a more structured conformation upon binding to RNA. We initiate simulations from a structured model of λ N in which the bound RNA is removed. The crystal structure corresponding to PDB IDs: 2OED (GB3), 1K61 (homeodomain), 1UBQ (ubiquitin), 7S6O (di-ubiquitin), 1XQ8 ( α -synuclein) and 1QFQ ( λN ) are shown in Figs. 2 (a), (b), (c) (d), (e) and (f) , respectively. Download figure Open in new tab Figure 2: Crystal structure of (a) GB3 (PDB ID:2OED), (b) homeodomain (PDB ID: 1K61), (c) ubiquitin (PDB ID: 1UBQ), (d) di-ubiquitin (PDB ID: 7S6O), (e) α -synuclein (PDB ID: 1XQ8) and (f) λ N (PDB ID:1QFQ) proteins. 2.2. All atom (AA) Molecular Dynamics (MD) simulations The GROMACS 50 2018.6 package 51 with the Amber99sb force field (ff) 52 is used for our AA simulation. The leapfrog algorithm is used to integrate the equations of motion. The TIP3P water model is used as the solvent. Periodic boundary conditions are applied in all three dimensions. The system is electrically neutralized by adding the required number of sodium (Na+) and chloride (Cl−) ions. The potential energy is minimized using the steepest descent algorithm. 53 Then AA MD simulation is performed at 300K temperature and 1 atmosphere pressure, maintaining an isothermal-isobaric (NPT) ensemble. We use the Berendsen thermostat 54 to maintain temperature and the Parrinello-Rahman barostat 55 to maintain constant pressure. The Lennard-Jones (LJ) and short-range electrostatic interactions are terminated at 10 Å. We use the Particle-Mesh Ewald (PME) 56 method to compute the long-range electrostatic interactions. LINCS 57 constraints are applied to all bonds involving hydrogen atoms. We use 2 fs time step for integration. The equilibration of the system is confirmed by the saturation of the root mean square deviation (RMSD) with time. We consider the equilibrated part of the trajectory for further analysis. Additionally, we simulate GB3 using the CHARMM27 force field (ff) 58 while keeping all other simulation conditions identical. 2.3. Coarse-grained (CG) Monte Carlo (MC) simulations The model system is simulated using the Metropolis Monte Carlo (MC) 42 algorithm in the canonical (NVT) ensemble. The model parameters are extracted from the AA trajectories. The system is maintained at a temperature k B T = 1, where k B is Boltzmann’s constant and T is the absolute temperature. A total of 100,000 MC steps are performed and the equilibration is monitored by observing the potential energy of the system. Post-equilibration trajectories are used to compute various observables, including solvent distribution around bead particles and secondary structure characteristics. To improve statistical reliability, multiple (5) independent simulations with identical initial configurations are conducted and the results are averaged. We use the diameter of the solvent bead, σ s (= 0.25 nm), as the unit of length and the thermal energy, k B T (= 2.5 kJ/mol), as the unit of energy in our simulations. We fix the solvent bead density at ~ 1 gm/cm 3 . 3. Analysis We calculate the following quantities over equilibrated trajectories of the AA MD and CG MC simulations. 3.1. Solvent distribution function We compute the solvent distribution function, 59 ρ ( r ) to understand the arrangement of solvent beads around solvophilic and solvophobic beads of the polymer. The ρ ( r ) is computed using the following function: Here, ⟨ ρ N ⟩ represents averaged number density of solvent beads around polymer bead, i index over solvent bead, j index over reference polymer bead, N is the total number of solvent beads, N 1 is the total number of solvophilic or solvophobic polymer beads, r ij is the distance between solvent bead i and polymer bead j, 4 π r 2 is the surface area of a spherical shell of radius r and δ ( r ij − r ) is the Dirac delta function. 3.2. Dihedral angle We compute the phi ( ϕ ) 60 and psi ( ψ ) 60 dihedral angles per residue per frame using our in-house program. 3.3. Structural persistence We calculate the structural persistence ( S P ) 14 , 61 parameter per residue using the formula: Where, N denotes the total number of frames. Δ ϕ i and Δ ψ i are the absolute values of the changes in dihedral angles ϕ and ψ of the residue in the frame i from the reference frame, Δ ϕ max and Δ ψ max are the maximum alterations possible in the Ramachandran diagram. 62 S P = 1 indicates no conformational change, whereas low S P represents greater deviation from the reference structure. 3.4. Ramachandran plot We generate the Ramachandran plot (RC) 62 using the averaged ϕ and ψ dihedral angles of residues over the equilibrated trajectory. 4. Results We test our coarse-grained (CG) model on a well-structured protein, GB3 and then apply this model to other well-structured proteins, like homeodomain, ubiquitin and di-ubiquitin. We also apply our methodology to intrinsically disordered proteins (IDPs), α -synuclein ( α S) and λ N. The resulting conformational ensembles are compared with available experimental crystal structures and all-atom (AA) simulation data. 4.1 Case of structured proteins We perform 1 µ s long AA simulation for GB3 starting from the crystal structure as the initial data. The RMSD over the course of the AA simulations is shown in SI Fig. S1 (a). The data show equilibration of the system. An equilibrium snapshot from AA simulations of the GB3 protein is shown in SI Fig. S1 (b). We take the bottom-up approach, namely using the AA data to build up the CG model. We use these CG parameters to perform MC simulations on various proteins. 4.1.A Building of the CG model parameters The COM of a residue in the crystal structure represents the center of the CG model bead. We treat the oxygen atom of the water as a solvent bead. We take polymer beads of two types: solvophilic and solvophobic, which repels the model solvent beads mimicking the hydrophilic and hydrophobic residues respectively. The beads representing the charged residues are taken to be solvophilic and the rest of the residues are taken to be solvophobic at a given pH. 63 We use the equilibrated AA trajectory of GB3 to compute different parameters for the CG model as detailed below. I. Bead size We compute the radius of gyration, R g of all the protein residues based on the distances of their heavy atoms from the COM. Construct the probability distribution of the radius of gyration, H ( R g ) considering all the residues over the equilibrium conformations. We assign free energy corresponding to the distribution, F r = − RT ln H ( Rg ), as shown in Fig. 3 (a) , where, R is the ideal gas constant and T is the temperature. This profile has minima which are the stable values of the radius of gyration in equilibrium. We assign the bead radius ( σ b / 2 ~ 0.2 nm), corresponding to the mean of the positions of the minima. Download figure Open in new tab Figure 3: (a) Free energy profile, − RT lnH ( R g ) vs radius of gyration of polymer bead, R g , (b) Free energy profile, − RT lnH ( d ) vs distance, d between two consecutive polymer beads and (c) Free energy profile, − RT lnH ( θ ) vs angle, θ between three consecutive polymer beads. II. Bonded interactions Stretching potential: We compute the distribution H ( d ) of the bond distance, d is the distance between center of mass of two consecutive amino acid residues. The stretching free energy, F s = − RT ln H ( d ), shown in Fig. 3 (b) . The harmonic bond stretching potential between two polymer beads is modeled as: , where r is the bond distance between two beads and s is the equilibrium bond distance given by the minimum of F s , shown in Fig. 3 (b) . The stretching force constant k s is calculated from the second-order derivative around the minimum of the stretching energy. Here, we find s = 0.5 nm and k s = 2292 kJ/mol/nm 2 . Bending potential: Fig. 3 (c) shows the bending free energy, F b = − RT ln H ( θ ), where H ( θ ) represents the distribution of the bond-angle θ between COM of three consecutive residues in AA data. The bending potential between three consecutive polymer beads is represented by a harmonic potential: , The equilibrium bond angle, θ 0 , is set to 1.4 rad, where F b has a minimum. The bending force constant, k a (87 kJ/mol/rad 2 ), is calculated from the second-order derivative around the minimum of − RT ln H ( θ ) vs θ in Fig. 3 (c) . III. Dihedral interactions The backbone dihedral angles ϕ and ψ are calculated from the atomic coordinates of the residues. We compute the joint probability distribution of the dihedral angles P (Γ i , Γ j ). Here, Γ i (= ϕ i , ψ i ) stands for backbone dihedral angles of the i th residue and Γ j (= ϕ j , ψ j ) are those for the j th residue. We define the dihedral interaction profile, F dih (Γ i , Γ j ) = − k B T ϵ ij lnP (Γ i , Γ j ). Here ϵ ij = 1 for i = j (intra-residual case) and the COM distance between i and j th residue ( i ≠ j , inter-residual case) is within a cut-off; otherwise, ϵ ij = 0. Let us first consider intra-residual distributions ( i = j ). We group the data into two distinct free energy profiles: one, considering the backbone dihedral angles of all the hydrophilic amino acid residues and the other, considering all the hydrophobic residues, shown in a two-dimensional grid, shown in Figs. 4 (a) and (b) , respectively. This classification takes care of the chemical properties of the side chains at the simplest level, although the side chains are not explicitly considered in our model. Download figure Open in new tab Figure 4: Free energy landscape for intra residual dihedral coupling, considering all (a) hydrophilic residues and (b) hydrophobic residues. A similar method is adopted to calculate the free energy profile for inter-residual ( i ≠ j ) dihedral coupling. The dihedral angles of the corresponding residues are subsequently classified according to the hydrophilic-hydrophilic, hydrophobic-hydrophobic and hydrophobic-hydrophilic residue pairs. For each case, all possible dihedral angle combinations are considered and the negative logarithm of each joint probability distribution yields the corresponding free energy profile. The free energy landscape (FEL) plots for various combinations of dihedral angles of two different residues are shown in Figs. 5 (a)-(l) . Download figure Open in new tab Figure 5: Free energy landscape for inter residual dihedral coupling considering (a-d) hydrophilic-hydrophilic residues, (e-h) hydrophilic-hydrophobic/hydrophobic-hydrophilic residues and (i-l) hydrophobic-hydrophobic residues. We, thus, have the dihedral angle interaction energy landscapes classified into solvophobic and solvophilic beads corresponding to hydrophobic and hydrophilic residues respectively. Moreover, we consider ( ϵ ij = 1) the inter-residue dihedral interaction energy only when the COM of the two beads are within a cut-off 7Å. Otherwise, we ignore the dihedral interaction energy ( ϵ ij = 0). The distinction between the free energy profiles between the intra and inter-residue dihedral angles within a cut-off distance implicitly take care of the interaction between the backbone and the dihedral angles. IV. Bead-solvent interaction We compute the distribution of water oxygen atoms ρ ( r ), where r is the distance of oxygen and the COM of the residue. ρ ( r ) data are averaged over conformations considering all the hydrophilic and hydrophobic residues separately. The data − RT ln( ρ ( r )) versus r for these two cases are shown in Figs. 6 (a) and (b) , respectively. Fig. 6 (a) shows the presence of a minimum at shorter distance, while no such minimum is observed in Fig. 6(b) . Download figure Open in new tab Figure 6: (a) Solvent density derived free energy profile, −RTln ρ ( r ) vs distance r , for water oxygen around hydrophilic residues, (b) Solvent density derived free energy profile, − RT lnρ ( r ) vs distance r , for water oxygen around hydrophobic residues, inset: Logarithmic free energy, ln F = − RT ln ρ ( r ), as a function of r for oxygen in water around hydrophobic residues and (c) Free energy profile, − RT lng ( r ), derived from the oxygen–oxygen radial distribution function, g ( r ) of water as a function of intermolecular distance, r. All the quantities are calculated from AA MD simulations data of GB3 protein. The minimum corresponding to hydrophilic residues is described via the Lennard-Jones (LJ) 12-6 potential, defining the solvophilic polymer bead-solvent bead interaction: . Here, r is the distance between the polymer bead and the solvent bead. The interaction parameter ϵ Hl is equal to 1.25 kJ/mol, which is the minimum of the energy, F Hl = − RT ln ρ ( r ) ( Fig. 6 (a) ). The hydrophobic beads interact with the solvent via the potential: V Hb-Solv (r) = 13.0 exp(−2.4 r ) (inset of Fig. 6 (b) ), obtained by fitting the semilog plot of the free energy, F Hb = − RT ln ρ ( r ) vs. r . V. Solvent-solvent bead interactions The oxygen-oxygen radial distribution function, g(r) 59 is calculated and the F Solv − Solv = − RT ln( g ( r )) versus r plot is shown in Fig. 6(c) . Here, R is the distance between the center of two oxygen atoms. For solvent-solvent beads, only the non-bonded Lennard-Jones (LJ) 12-6 potential is considered: . The interaction parameter ϵ SS , equal to 2.5 kJ/mol, corresponds to the minimum value of the energy, F Solv − Solv ( Fig. 6(c) ). VI. Non-bonded bead-bead interaction Repulsive interaction: , where ϵ − is the interaction parameter, σ b is the diameter of the polymer bead and r ij is the distance between two polymer beads i and j . This is to ensure that the beads cannot overlap. We set the ϵ bb value to 2.5 kJ/mol, given by the minimum of the F r ( Fig. 3 (a) ). Screened Coulomb interactions: , where q and q are the net charges of the residues i and j , respectively, as calculated using the pdb2pqr serve 64 at a given pH value. We take pH = 7 in our calculation. ϵ 0 is the permittivity of vacuum and D is the dielectric constant of water (set equal to 80). Thus, medium for the electrostatic interaction is described as a continuum. We set the Debye screening length, λ equal to 0.03 nm computed from the salt concentration and charge of the residue. Here, e is the elementary charge, z i (= 1) the valency of i type of ion present in the solution, c i the ionic concentration (0.1 mol L −1 ), and N A the Avogadro’s number. 4.1.B Monte Carlo for the CG model The CG model potential that considers different interaction terms in the system takes the following form: Here, r ij is the distance between two polymer beads having coordinates and denotes the bond angle between three consecutive polymer beads i, j , and k . is the coordinate of the α th solvent bead. The prime over the bonding term indicates that only two consecutive residues COM and the double prime over the bending energy term indicates that three consecutive residues are to be taken. for solvophilic beads and for solvophobic bead at and solvent bead is the distance between two solvent beads are having coordinates and . The bead-dihedral interactions are implicitly accounted for by bead-bead distance cut-off where the dihedral interactions are considered for beads within the cut-off distances. For the inter-bead dihedral contribution, we consider only those beads that are within 7.0 Å of each bead in every frame. The bead-solvent and solvent-solvent interactions are truncated at a cutoff distance of 6.25 Å . We take all the parameters in the model potential as determined from the GB3 AA trajectory, given in Table 1 . The initial structure is taken from the crystal structure of GB3. Each bead is located at the center of mass of a residue determined from the heavy atoms in a given residue. The beads are classified as either solvophilic or solvophobic, as given in SI Table S1, depending on hydrophilic or hydrophobic residue obtained from water oxygen distributions. The initial values of the dihedral angles of different residues are taken from the crystal structure and the solvent molecules are the oxygen atoms of water molecules in an equilibrated AA configuration. The coarse-grained (CG) polymer and solvent beads are placed in a cubic simulation box of length L=5.6 nm, with the periodic boundary conditions applied in all directions (x, y, and z). We use 56 polymer beads and 5,444 solvent beads in our simulation. View this table: View inline View popup Download powerpoint Table 1: Table of CG parameters. During a MC move, a bead, either polymer or solvent, is selected at random. All three position coordinates along with two dihedral angles are given a random movement to generate a new configuration of the selected polymer bead. For a solvent bead, only three position coordinates are updated. The energy changes due to the change in the degrees of freedom of the bead are computed based on various interactions in the CG potential, V CG . The contribution due to changes in the dihedral angles, V dih (Γ i , Γ j ) is computed as follows. The intra-bead (i=j) interaction between the dihedral of the beads is calculated by linear interpolation from the two-dimensional free energy profile, F grids in Fig. 4 . We choose the type of grid depending on if the bead is solvophilic or solvophobic. The inter-bead (i≠ j) dihedral coupling depends on the solvophilic or solvophobic character of the pair of beads. We use the V dih (Γ i , Γ j ) parameters shown in Fig. 5 to interpolate. We allow the system to equilibrate, which is monitored via the total energy (SI Fig. S2) and calculate the structural quantities over the equilibrated trajectory. We calculate the solvent bead distributions ρ ( r ) at a distance r from the bead center around solvophilic and solvophobic beads ( Fig. 7 (a) ) in the scaled unit. Solvophilic beads form strong interactions with solvent beads, leading to the formation of a well-defined first solvation shell. As a result, the local solvent density increases near the reference bead, producing a pronounced first peak in solvent distribution. In contrast, solvophobic beads do not interact strongly with the solvent and therefore the solvent distribution near the bead is not significantly enhanced; consequently, the first peak is weak or absent. In the inset, the corresponding AA data are shown. For the AA results, the x-axis represents the distance between the COM of the given residue (hydrophilic or hydrophobic) and oxygen atom of water particles in nm, while the y-axis represents the normalized average solvent number density at distance r , which is also dimensionless. Similar to the CG solvophilic beads, hydrophilic residues in the AA data form strong interactions with solvent molecules and stabilize structured hydration shells, leading to a clear first peak in the solvent distribution. On the other hand, hydrophobic residues interact weakly with solvent molecules, so the first peak in the solvent density is weak or absent. Overall, the CG results reproduce the qualitative trend observed in the AA data. We further compute structural persistence, S P per residue in equilibrated CG MC and AA MD trajectories, as shown in Fig. 7 (b). The CG and AA simulation data for S P exhibit good agreement. Download figure Open in new tab Figure 7: (a) Solvent distribution ρ ( r ) around solvophilic (salmon, solid line) and solvophobic (turquoise, dashed line) beads in the equilibrated CG MC trajectory of GB3. The inset shows the solvent distribution around hydrophilic (salmon, solid line) and hydrophobic (turquoise, dashed line) residues in the equilibrated AA MD trajectory and (b) Structural persistence (S P ) per residue in equilibrated CG simulations (salmon, solid line) and AA MD simulations (turquoise, dashed line). We examine how far the structural aspects are captured in our CG simulations. To this end, we compare the secondary structure preferences of each residue from the CG MC simulations. We assign the secondary structure elements to each bead, namely, helix (H), β -sheet (S) and unstructured (U) depending on the values ϕ and ψ for each conformation. 65 The unstructured (U) category includes all conformations that do not fall into the helix (H) or β -sheet (S). We assign a given bead the preferred structural element, the one that occurs the maximum number of times over the equilibrium trajectory. A similar strategy is followed for the AA data and the secondary structural elements are assigned directly from the dihedral angles of the backbone. SI Table S2 presents the structural preferences of each residue in CG MC simulation, along with their corresponding secondary structures in the crystal structure and AA MD simulations. In Fig. 8(a) and (b) , we show the probability distributions of the ϕ and ψ dihedral angles for LYS10 in equilibrated trajectories. The distributions from AA MD and CG MC simulations show good agreement. We find that the secondary structural elements show ~80% agreement with the residues observed in the crystal structure and AA data. Download figure Open in new tab Figure 8: Probability distributions of the (a) ϕ and (b) ψ dihedral angles of LYS10 in the GB3 protein, calculated from equilibrated trajectories. CG MC simulations are shown in salmon, while AA MD simulations are shown in turquoise. (c) Comparison of Ramachandran (RC) plots of the GB3 protein for the crystal structure (green), the final structure from all-atom molecular dynamics (AA MD) simulations (red), and the final structure from coarse-grained Monte Carlo (CG MC) simulations (black). Download figure Open in new tab Figure 9: (a) Bar plot showing residue counts in % in helix, sheet, and unstructured regions for five cases—AA simulation (AA), our structure-based CG force field (CG ff), SIRAH, Martini, and the crystal structure. (b) Secondary-structure similarity matrix, where off-diagonal elements quantify agreement between methods. We further analyze the Ramachandran (RC) plots to compare the final structures obtained from the CG simulations with those from the AA simulations and the crystal structure of the GB3 protein. The RC plots, shown in Fig. 8 (c) , confirms excellent agreement. We show in Table 2 comparison of percentages of residues in a given structural element over the CG and AA trajectories along with those in the crystal structure. the helix (H) elements are in excellent agreement to the AA and crystal structure data. Our model slightly over-estimates the β -sheet (S) compared to the crystal structure, while the AA data underestimates the β -sheet (S) structure. The unstructured (U) region in the CG model agrees better with the crystal structure compared to the AA data. Overall, our model reproduces the secondary structural elements as accurately as in the AA data. View this table: View inline View popup Download powerpoint Table 2: Comparison of secondary structure preferences of proteins. AA denotes trajectory of all atom molecular dynamics simulations and CG denotes conformations based on coarse-grained Monte Carlo simulations. ‘H’ corresponds Helix, ‘S’ corresponds sheet and ‘U’ corresponds to other than element helix or sheet i.e. loop/coil/turn/bend region of the protein. We carry out controlled studies to check the robustness of our model. It is unclear whether the good agreement between the CG and crystal structures arises from any bias introduced by our initial choices. To investigate this, for GB3, we add to the crystal dihedral angles random number within a range of − π/ 10 to + π/ 10. We perform the CG simulations starting from the random structure. Here also, the resulting residue-wise secondary structure matches well with the crystal structure (80.3%) (SI Table S3). As shown in Table 2 , our CG model reproduces the percentages of different structural elements in very good agreement to the AA and the crystal structure data. We examine whether using the AA trajectory introduces any force field bias in the CG parameters. To check this, we perform CG simulations on GB3 model system using interaction parameters derived from widely used AA force field, CHARMM and compare the resulting structural features with crystal data and those obtained from corresponding AA MD simulations. The data are shown in SI Tables S4. The structural agreement observed in all cases is consistent with our previous findings based on the AMBER force field, showing > =75% similarity with the crystal and AA data (SI Tables S4), along with similar percentages of different secondary structural elements ( Table 2 ). 4.1.C Transferability of the CG model parameters to other structured proteins Next, we check how far the structural data for other well-structured proteins can be reproduced by using the CG parameters derived from GB3 AA trajectory. We consider three different well structured proteins, namely, homeodomain, ubiquitin and di-ubiquitin. For each of them we generate structural data from AA trajectories and CG simulations with CG parameters as those of the GB3 protein starting from crystal structures and using the AMBER force-field. The charged residues in these proteins are taken to be solvophilic, while the other beads are classified to be solvophobic. SI Table S5 shows the secondary structural elements of each residue derived from the crystal structure, AA simulations and CG MC simulations for homeodomain. We observe good agreement: The CG MC simulation reproduces the crystal structural data with 90% accuracy (SI Table S5). The residue-wise secondary structural elements from CG MC simulations agree well with the AA MD results, showing an agreement of 81% (SI Table S5). Our CG model captures the secondary structural elements as good as the AA data and are quite comparable to the crystal structure data ( Table 2 ). In case of ubiquitin (PDB ID : IUBQ), we find that both the AA and CG data overestimates the structured portion of the protein, while underestimating the unstructured parts. SI Table S6 shows detailed comparison: The secondary structural elements show 67% matching with crsytal structure. Although this agreement is poorer compared to other proteins, the mismatched residues mostly correspond to the unstructured (U) conformation in the crystal structure. The residues adopting structured conformations (H or S) determined from CG simulations agree well with the crystal structure. Compared to the AA MD simulation data (also in SI Table S6), the agreement improves significantly to 92%. Table 2 lists the agreements secondary structural element wise: Here also, the structure elements in particular H conformation is overestimated in both AA and CG models, while the U elements are underestimated compared to the crystal structure. We also study di-ubiquitin in which two ubiquitin molecules are connected by a covalent isopeptide bond (PDB ID: 7S6O), where the C-terminal glycine (Gly76) of the first ubiquitin is linked to lysine (Lys48) of the second ubiquitin. We compare residue-wise secondary structure preferences obtained from the CG simulations with those derived from the crystal structure and AA simulations (SI Tables S7 (a) & (b)). The CG results show good agreement with the crystal structure, approximately 70% of residues exhibiting matching secondary structure preferences. Most of the mismatched residues are having a U conformation in the crystal structure. Compared to the AA MD data, approximately 74% of residues show consistent secondary structure preferences, while the remaining residues predominantly adopt S or U conformations ( Table 2 ). Table 2 shows a comparative analysis of percentage of secondary structural elements of di-ubiquitin. The data show that the agreement of the secondary structural elements is quite good for all the cases. 4.2. Case of intrinsically disordered protein (IDP) Next, we consider the cases of couple of IDPs: α S and λ N. The initial structure of α S is taken from the crystal structure deposited under PDB ID: 1XQ8. We compute all bonded and non-bonded interaction parameters from the AA MD trajectory of α S and find that the bond force constant ( k b = 746 kJ mol −1 nm −2 , as detailed in the SI) and the angle force constant ( k a = 30 kJ mol −1 rad −2 , as detailed in the SI) are smaller than those obtained for GB3, while the remaining parameters are comparable. These reduced elastic constants reflect the increased conformational flexibility inherent to the intrinsically disordered nature of the α S protein. We perform CG simulations of α S using these force constants and compare residue-wise secondary structure preferences with those obtained from the AA MD simulations data, as shown in SI Tables S8 (a) & (b). Since the IDPs lack well defined crystal structure, we compare the structural elements from our CG simulations to those from AA data only. Comparison with the AA MD data shows that approximately 67.4% of residues exhibit consistent secondary structure preferences (SI Tables S8 (a) & and (b)). In this case, the unmatched residues in the CG simulations predominantly sample U conformations, whereas they adopt H conformations in the AA simulations. Table 2 shows the overall percentages of the secondary structure elements in the CG and AA data. The U elements are in vast majority in CG structures than the AA data. The lack of structure is consistent with the intrinsically disordered nature of the protein which is better reproduced in the CG simulations than that in the AA simulations. We adopt the parameter set derived from α S for the CG simulations of model λ N. Here also we observe that in both CG and AA simulations, the residues preferentially adopt U conformations ( Table 2 and SI Table S9) as expected for IDPs. 5. Discussions An explicit comparison of our model to other widely used CG models is worthwhile here. We perform to this end simulations of the GB3 protein using both the SIRAH and Martini CG force fields in explicit solvent. The details of these simulation methods are provided in the SI. The final GB3 structures after 1 µ s simulations using SIRAH and Martini CG force fields are shown in SI Fig.S3 (a) and (b), respectively, with residues colored by secondary structure (helix, sheet, unstructured). Next, we count the secondary-structural elements based on final structure for four cases: AA simulation, our CG force field, SIRAH, Martini, as well as experimental crystal structure—to provide a clear basis for comparison (see Fig.9a ). The CG force field SIRAH exhibits a bias toward disorder, underestimating helices and sheets. In contrast, the Martini force field shows the opposite tendency i.e. over-stabilizing ordered secondary structure and significantly under predicting unstructured structures compared to experiment. on the other hand, our CG model is able to show results closer to both AA and crystal structure data. We further construct similarity matrix by performing pair wise residue-by-residue comparison of secondary structure assignments (helix, sheet or unstructured) between all five systems ( Fig.9(b) ). For each pair of system, the percentage similarity has been calculated as the number of residues with identical secondary structure assignments divided by the total number of residues, multiplied by 100. It gives a symmetric matrix with 100% self-similarity on the diagonal and off-diagonal elements quantify the agreement between different methodologies. Our CG model shows strongest overall performance, achieving the highest pairwise similarity with the AA simulation (81.8%), indicating it most faithfully reproduce the structural details captured by the AA reference. Our CG model also maintain a good agreement with experimental crystal structure as well (74.5%). On the other hand, Martini shows relatively good agreement with the crystal structure (80.0%) but poor agreement with AA (63.4%) and SIRAH (65.5%). Conversely, SIRAH’s strongest agreement is with AA (76.4%), but it shows lower similarity to crystal (70.9%) and the lowest agreement with Martini (65.5%). The matrix thus quantitatively confirm that our CG force field simultaneously preserves AA-level structural detail and maintains good agreement with the experimental fold. It therefore provides the most balanced and reliable residue-level secondary structure representation among the tested CG models. Earlier CG models for DNA–protein systems 66 , 67 and DPG-MC approach to generate polypeptides 39 and protein structures 40 heavily rely on data bank structures. Such approaches have inherent limitations: The data bank do not generate information on interaction between water and the bio-molecules, although these interactions play primary roles in adopting bio-molecular structure. The coarse-graining requires statistical averages over a wide range of conformations to accurately derive effective interaction potentials. It is difficult to obtain adequate averaging using limited number of available data bank structures. Moreover, the PDB structures lack consistent ensemble information and may be influenced by crystallization artifacts. Our current CG model surpasses these limitations by deriving the CG parameters from AA trajectories. This ensures that the solvent interactions are taken care of, the thermodynamic ensemble is well defined and there are sufficient conformation to extract meaningful CG parameters. We find, moreover, that the CG parameters obtained from different AA models can reproduce the structure with comparable accuracy. Our computational model offers several advantages compared to other CG models reported in the literature: (i) It is a single-site model. (ii) In each MC step, both Cartesian coordinates and backbone dihedral angles are updated simultaneously, avoiding separate Cartesian and dihedral-space sampling schemes and resulting in improved computational efficiency. (iii) It employs an explicit-solvent representation and distinguishes the beads based on solvent interaction. (iv) The backbone hydrogen bonding interactions are incorporated implicitly through the CG interaction profiles derived from atomistic simulation trajectories. The CG force field parameters derived from GB3 are transferable to other structured proteins. However, when we use the stretching and bending force constants derived from the GB3 protein, we observe lack of equilibration in both IDP proteins. This suggests that the backbone flexibilities play key role in adopting the equilibrium structure. On the other hand, the CG parameters for one IDP work well for another IDP. The transferability of the system arises from two main reasons: (1) The chemical environment of the backbone formed from peptide bonds is similar for different proteins with comparable elastic constants. (2) A broad categorization of solvent and dihedral interactions based on hydrophilic and hydrophobic residues is sufficient to describe the protein secondary structural elements. However, it must be pointed out that we have tested transferability only for restricted number of cases. A general proof of transferability is quite difficult. The main importance of our study is as follows: Given AA trajectory for a protein, we show a method to carry out systematic coarse graining while retaining the structural information. Based on CG model the phase behaviour and dynamics of an aqueous suspension of proteins can be studied which is beyond the scope of AA simulations. Our model can be improved in several important ways: (1) The electrostatic interaction effects are included as a constant with dielectric constant ϵ , considering the medium as a continuum. On the other hand, we use a monoatomic water model, where each oxygen atom of a water molecule is treated as a solvent bead to calculate the solvent-solvent and bead-solvent interactions. It is more important to treat all the interactions on equal footing. In particular, taking the dipole moment of water molecules explicitly into account would lead to a more realistic solvent model instead of using dielectric continuum.(2) It may be possible to include model terms to describe hydrogen bonding between water and the side chains to improve performance of the model to predict function. (3) Inclusion of side chain information is needed to describe ligand binding interactions through side chains. 6. Conclusion To summarize, we show that the structure of the protein can be captured by introducing appropriate CG model potential energy with backbone elasticity, dihedral interaction and distinction of solvent interaction with solvopobic and solvophilic beads. The CG force-field parameters are derived from AA simulation trajectories. This model would enable to perform much faster simulations for studying phenomena involving multiple protein molecules, such as protein aggregation. Our CG model needs to be improved to incorporate side chain dihedrals and water dipoles moment and can be extended study the in bio-molecular systems in terms of the CG variables. Similar CG model may extended to other bio-molecular systems. Data Availibility We have used GROMACS 2018 MD simulation package to perform all the atomistic simulations. The software can be found at https://manual.gromacs.org/documentation/respectively . The coarse grained monte carlo simulation has been using our in house code available at https://github.com/snbsoftmatter/mc-protein-cg . Author information Authors and Affiliations Kanika Kole -Department of Physics of Complex Systems, S. N. Bose National Centre for Basic Sciences, Block-JD, Sector-III, Salt Lake, Kolkata-700106, India. Abhik Ghosh Moulick - Department of Physics of Complex Systems, S. N. Bose National Centre for Basic Sciences, Block-JD, Sector-III, Salt Lake, Kolkata-700106, India. Jaydeb Chakrabarti - Department of Chemical and Biological Sciences, S. N. Bose National Centre for Basic Sciences, Block-JD, Sector-III, Salt Lake, Kolkata-700106, India. Author contributions KK and AGM: Conceptualization, data curation, formal analysis, methodology, validation, visualization, software, investigation, writing – original draft, writing – review & editing. J. Chakrabarti: conceptualization, methodology, validation, project administration, resources, supervision, writing – review & editing. Ethics declarations Competing interests The authors declare no competing interests. Acknowledgement The authors thank the Technical Research Centre, S. N. Bose National Centre for Basic Sciences for computational support. K.K. and A.G.M. thank Anirban Paul for helpful discussions. K.K. thanks the University Grants Commission (UGC), India [F. No. 16-6(DEC. 2018)/2019(NET/CSIR)], and A.G.M. thanks DST, India, for an INSPIRE fellowship (IF170961) as financial support. JC acknowledges financial support from the CSIR Emeritus Scientist scheme. Footnotes E-mail: kanikakole0094{at}gmail.com , ghoshabhik.physics{at}gmail.com (c) In the text 1. We have modified the Introduction section (Pages 2 to 6) by adding a brief discussion of earlier approaches, the need for a new approach, and the motivation and significance of applying the current coarse-grained modeling to IDPs. 2. We have added two new systems, di ubiquitin and α-synuclein, in the System Preparation subsection (Page 7). 3. We have added detailed simulation conditions for the CHARMM force field in the All-atom Simulation subsection (Page 9). 4. We have removed the Dihedral Principal Component Analysis subsection from the Analysis section. 5. We have removed the GROMOS force field simulation results. 6. We have incorporated the simulation results of the multi-domain protein di-ubiquitin (Page 28). 7. We have added the simulation results of an IDP, α synuclein (Page 29). 8. We have added a comparison of the results of our coarse-grained (CG) model with other models, Martini and SIRAH, in the Discussion section (Pages 30 31). 9. We have added the advantages of our CG model compared with other CG models in the Discussion section (Page 32). 10. We have modified the Conclusion section. (d) In the figure 1. We have added two new figures in Fig. 2 (Page 8). 2. We have split Fig. 3 into Figs. 3 (Page 14) and 6 (Page 18). 3. We have modified Fig. 7 (Page 23) (Fig. 6 in the previous version). 4. We have modified Fig. 8 (Page 25), (Fig. 7 in the previous version). 5. We have added Fig. 9 (Page 31). References (1). ↵ Saunders , M. G. ; Voth , G. A. Coarse-graining methods for computational biology . Annual review of biophysics 2013 , 42 , 73 – 93 . OpenUrl PubMed (2). ↵ Souza , P. C. T. ; Borges-Araújo , L. ; Brasnett , C. ; Moreira , R. A. ; Grünewald , F. ; Park , P. ; Thallmair , S. GōMartini 3: From Large Conformational Changes in Proteins to Environmental Bias Corrections . Nature Communications 2025 , 16 , 4051 . OpenUrl PubMed (3). Klein , M. L. ; Shinoda , W. Large-Scale Molecular Dynamics Simulations of Self-Assembling Systems . Science 2008 , 321 , 798 – 800 . OpenUrl Abstract / FREE Full Text (4). Voth , G. A. Coarse-Graining of Condensed Phase and Biomolecular Systems ; CRC Press : Boca Raton , 2008 . (5). Noid , W. G. Perspective: Coarse-Grained Models for Biomolecular Systems . The Journal of Chemical Physics 2013 , 139 , 090901 . OpenUrl CrossRef PubMed (6). Brini , E. ; Algaer , E. A. ; Ganguly , P. ; Li , C. ; Rodriguez-Ropero , F. ; van der Vegt , N. F. A. Systematic Coarse-Graining Methods for Soft Matter Simulations – A Review . Soft Matter 2013 , 9 , 2108 – 2119 . OpenUrl (7). ↵ Ingólfsson , H. I. ; López , C. A. ; Uusitalo , J. J. ; de Jong , D. H. ; Gopal , S. M. ; Periole , X. ; Marrink , S. J. The Power of Coarse Graining in Biomolecular Simulations . Wiley Interdisciplinary Reviews: Computational Molecular Science 2014 , 4 , 225 – 248 . OpenUrl PubMed (8). ↵ De Jong , D. H. ; Singh , G. ; Bennett , W. D. ; Arnarez , C. ; Wassenaar , T. A. ; Schäfer , L. V. ; Periole , X. ; Tieleman , D. P. ; Marrink , S. J. Improved Parameters for the Martini Coarse-Grained Protein Force Field . Journal of Chemical Theory and Computation 2013 , 9 , 687 – 697 . OpenUrl (9). ↵ Alessandri , R. ; Souza , P. C. T. ; Thallmair , S. ; Melo , M. N. ; De Vries , A. H. ; Marrink , S. J. Pitfalls of the Martini Model . Journal of Chemical Theory and Computation 2019 , 15 , 5448 – 5460 . OpenUrl (10). ↵ Charron , N. E. et al. Navigating protein landscapes with a machine-learned transferable coarse-grained model . Nature Chemistry 2025 , 17 , 1284 – 1292 . OpenUrl PubMed (11). ↵ IUPAC Compendium of Chemical Terminology, 5th ed. (the “Gold Book”) . doi: 10.1351/goldbook.D01730 , 2025 ; Online version (2006–): “Dihedral angle”. OpenUrl CrossRef (12). ↵ Habchi , J. ; Tompa , P. ; Longhi , S. ; Uversky , V. N. Introducing protein intrinsic disorder . Chemical reviews 2014 , 114 , 6561 – 6588 . OpenUrl CrossRef PubMed Web of Science (13). ↵ Dunker , A. K. et al. Intrinsically Disordered Protein . Journal of Molecular Graphics and Modelling 2001 , 19 , 26 – 59 . OpenUrl CrossRef (14). ↵ Moulick , A. G. ; Chakrabarti , J. Conformational Fluctuations in the Molten Globule State of α-Lactalbumin . Physical Chemistry Chemical Physics 2022 , 24 , 21348 – 21357 . OpenUrl PubMed (15). ↵ Moulick , A. G. ; Chakrabarti , J. Fluctuation-Dominated Ligand Binding in Molten Globule Protein . Journal of Chemical Information and Modeling 2023 , 63 , 5583 – 5591 . OpenUrl PubMed (16). ↵ Asthana , S. ; Bhattacharyya , D. ; Kumari , S. ; Nayak , P. S. ; Saleem , M. ; Bhunia , A. ; Jha , S. Interaction with Zinc Oxide Nanoparticle Kinetically Traps α-Synuclein Fibrillation into Off-Pathway Non-Toxic Intermediates . International Journal of Biological Macromolecules 2020 , 150 , 68 – 79 . OpenUrl PubMed (17). ↵ Takada , S. Coarse-Grained Molecular Simulations of Large Biomolecules . Current Opinion in Structural Biology 2012 , 22 , 130 – 137 . OpenUrl CrossRef PubMed (18). ↵ Dans , P. D. ; Zeida , A. ; Machado , M. R. ; Pantano , S. A coarse grained model for atomicdetailed DNA simulations with explicit electrostatics . Journal of chemical theory and computation 2010 , 6 , 1711 – 1725 . OpenUrl (19). ↵ Marrink , S. J. ; Risselada , J. H. ; Yefimov , S. ; Tieleman , D. P. ; de Vries , A. H. The MARTINI Force Field: Coarse Grained Model for Biomolecular Simulations . Journal of Physical Chemistry B 2007 , 111 , 7812 – 7824 . OpenUrl CrossRef (20). ↵ Barrera , E. E. ; Machado , M. R. ; Pantano , S. Fat SIRAH: coarse-grained phospholipids to explore membrane–protein dynamics . Journal of Chemical Theory and Computation 2019 , 15 , 5674 – 5688 . OpenUrl (21). ↵ Lopez , C. A. ; Rzepiela , A. J. ; de Vries , A. H. ; Dijkhuizen , L. ; Hünenberger , P. H. ; Marrink , S. J. Martini Coarse-Grained Force Field: Extension to Carbohydrates . Journal of Chemical Theory and Computation 2009 , 5 , 3195 – 3210 . OpenUrl (22). ↵ Darré , L. ; Machado , M. R. ; Dans , P. D. ; Herrera , F. E. ; Pantano , S. Another coarse grain model for aqueous solvation: WAT FOUR? Journal of Chemical Theory and Computation 2010 , 6 , 3793 – 3807 . OpenUrl (23). ↵ Hadley , K. R. ; McCabe , C. Coarse-Grained Molecular Models of Water: A Review . Molecular Simulation 2012 , 38 , 671 – 681 . OpenUrl PubMed (24). ↵ Dill , K. A. Theory for the Folding and Stability of Globular Proteins . Biochemistry 1985 , 24 , 1501 – 1509 . OpenUrl CrossRef PubMed Web of Science (25). Lau , K.-F. ; Dill , K. A. A Lattice Statistical Mechanics Model of the Conformational and Sequence Spaces of Proteins . Macromolecules 1989 , 22 , 3986 – 3997 . OpenUrl CrossRef Web of Science (26). Dill , K. A. ; Fiebig , K. M. ; Chan , H. S. Cooperativity in Protein-Folding Kinetics . Proceedings of the National Academy of Sciences of the United States of America 1993 , 90 , 1942 – 1946 . OpenUrl Abstract / FREE Full Text (27). ↵ Farris , A. C. ; Shi , G. ; Wüst , T. ; Landau , D. P. The role of chain-stiffness in lattice protein models: A replica-exchange Wang–Landau study . J. Chem. Phys . 2018 , 149 , 125101 . OpenUrl PubMed (28). ↵ Kmiecik , S. ; Gront , D. ; Kolinski , M. ; Wieteska , L. ; Dawid , A. E. ; Kolinski , A. Coarse-Grained Protein Models and Their Applications . Chem. Rev . 2016 , 116 , 7898 – 7936 . OpenUrl CrossRef PubMed (29). Davtyan , A. ; Schafer , N. P. ; Zheng , W. ; Clementi , C. ; Wolynes , P. G. ; Papoian , G. A. AWSEM-MD: protein structure prediction using coarse-grained physical potentials and bioinformatically based local structure biasing . The Journal of Physical Chemistry B 2012 , 116 , 8494 – 8503 . OpenUrl PubMed (30). ↵ Mukherjee , A. ; Bagchi , B. Correlation between rate of folding, energy landscape, and topology in the folding of a model protein HP-36 . J. Chem. Phys . 2003 , 118 , 4733 – 4747 . OpenUrl (31). ↵ Darré , L. ; Machado , M. R. ; Brandner , A. F. ; González , H. C. ; Ferreira , S. ; Pantano , S. SIRAH: A Structurally Unbiased Coarse-Grained Force Field for Proteins with Aqueous Solvation and Long-Range Electrostatics . J. Chem. Theory Comput . 2015 , 11 , 723 – 739 . OpenUrl CrossRef PubMed (32). ↵ Ghosh Moulick , A. ; Patel , R. ; Onyema , A. ; Loverde , S. M. Unveiling Nucleosome Dynamics: A Comparative Study Using All-Atom and Coarse-Grained Simulations Enhanced by Principal Component Analysis . J. Chem. Phys . 2025 , 162 , 065101 . OpenUrl PubMed (33). ↵ Kharche , S. ; Yadav , M. ; Hande , V. ; Prakash , S. ; Sengupta , D. Improved Protein Dynamics and Hydration in the Martini3 Coarse-Grain Model . J. Chem. Inf. Model . 2024 , 64 , 837 – 850 . OpenUrl CrossRef PubMed (34). ↵ Monticelli , L. ; Kandasamy , S. K. ; Periole , X. ; Larson , R. G. ; Tieleman , D. P. ; Marrink , S. J. The MARTINI Coarse-Grained Force Field: Extension to Proteins . J. Chem. Theory Comput . 2008 , 4 , 819 – 834 . OpenUrl CrossRef PubMed Web of Science (35). ↵ Borges-Araújo , L. ; Patmanidis , I. ; Singh , A. P. ; Santos , L. H. ; Sieradzan , A. K. ; Vanni , S. ; Souza , P. C. T. Pragmatic Coarse-Graining of Proteins: Models and Applications . J. Chem. Theory Comput . 2023 , 19 , 7112 – 7135 . OpenUrl CrossRef PubMed (36). ↵ Kawamoto , S. ; Liu , H. ; Miyazaki , Y. ; Seo , S. ; Dixit , M. ; DeVane , R. ; MacDermaid , C. ; Fiorin , G. ; Klein , M. L. ; Shinoda , W. SPICA force field for proteins and peptides . Journal of Chemical Theory and Computation 2022 , 18 , 3204 – 3217 . OpenUrl (37). ↵ Liwo , A. ; Ołdziej , S. ; Pincus , M. R. ; Wawak , R. J. ; Rackovsky , S. ; Scheraga , H. A. A united-residue force field for off-lattice protein-structure simulations. I. Functional forms and parameters of long-range side-chain interaction potentials from protein crystal data . Journal of computational chemistry 1997 , 18 , 849 – 873 . OpenUrl CrossRef Web of Science (38). ↵ Wassenaar , T. A. ; Pluhackova , K. ; Bockmann , R. A. ; Marrink , S. J. ; Tieleman , D. P. Going backward: a flexible geometric approach to reverse transformation from coarse grained to atomistic models . Journal of chemical theory and computation 2014 , 10 , 676 – 690 . OpenUrl (39). ↵ Evans , J. S. ; Chan , S. I. ; Mathiowetz , A. M. ; Goddard , W. A. I. De Novo Prediction of Polypeptide Conformations Using Dihedral Probability Grid Monte Carlo Methodology . Protein Sci . 1995 , 4 , 1203 – 1216 . OpenUrl PubMed Web of Science (40). ↵ Mathiowetz , A. M. ; Goddard , W. A. I. Building Proteins from Cα Coordinates Using the Dihedral Probability Grid Monte Carlo Method . Protein Sci . 1995 , 4 , 1217 – 1232 . OpenUrl PubMed Web of Science (41). ↵ Liang , F. ; Wong , W. H. Evolutionary Monte Carlo for Protein Folding Simulations . J. Chem. Phys . 2001 , 115 , 3374 – 3380 . OpenUrl CrossRef (42). ↵ Theodorou , D. N. Progress and Outlook in Monte Carlo Simulations . Industrial & Engineering Chemistry Research 2010 , 49 , 3047 – 3058 . OpenUrl (43). ↵ Kociurzynski , R. ; Makshakova , O. N. ; Knecht , V. ; Romer , W. Multiscale Molecular Dynamics Studies Reveal Different Modes of Receptor Clustering by Gb3-Binding Lectins . J. Chem. Theory Comput . 2021 , 17 , 2488 – 2501 . OpenUrl PubMed (44). ↵ Ulmer , T. S. ; Ramirez , B. E. ; Delaglio , F. ; Bax , A. Evaluation of Backbone Proton Positions and Dynamics in a Small Protein by Liquid Crystal NMR Spectroscopy . J. Am. Chem. Soc . 2003 , 125 , 9179 – 9191 . OpenUrl CrossRef PubMed Web of Science (45). ↵ Aishima , J. ; Gitti , R. K. ; Noah , J. E. ; Gan , H. H. ; Schlick , T. ; Wolberger , C. A Hoogsteen Base Pair Embedded in Undistorted B-DNA . Nucleic Acids Res . 2002 , 30 , 5244 – 5252 . OpenUrl CrossRef PubMed Web of Science (46). ↵ Vijay-Kumar , S. ; Bugg , C. E. ; Cook , W. J. Structure of Ubiquitin Refined at 1.8 Å Resolution . J. Mol. Biol . 1987 , 194 , 531 – 544 . OpenUrl CrossRef PubMed Web of Science (47). ↵ Wydorski , P. M. ; Osipiuk , J. ; Lanham , B. T. ; Tesar , C. ; Endres , M. ; Engle , E. ; Jedrzejczak , R. ; Mullapudi , V. ; Michalska , K. ; Fidelis , K. ; Fushman , D. ; Joachimiak , A. ; Joachimiak , L. A. Dual Domain Recognition Determines SARS-CoV-2 PLpro Selectivity for Human ISG15 and K48-Linked Di-Ubiquitin . Nature Communications 2023 , 14 , 1 – 16 . OpenUrl PubMed (48). ↵ Ulmer , T. S. ; Bax , A. ; Cole , N. B. ; Nussbaum , R. L. Structure and Dynamics of Micelle-Bound Human alpha-Synuclein . Journal of Biological Chemistry 2005 , 280 , 9595 – 9603 . OpenUrl Abstract / FREE Full Text (49). ↵ Schärpf , M. ; Sticht , H. ; Schweimer , K. ; Boehm , M. ; Hoffmann , S. ; Rösch , P. Antitermination in Bacteriophage λ: The Structure of the N36 Peptide-boxB RNA Complex . Eur. J. Biochem . 2000 , 267 , 2397 – 2408 . OpenUrl CrossRef PubMed Web of Science (50). ↵ Abraham , M. J. ; Murtola , T. ; Schulz , R. ; Pall , S. ; Smith , J. C. ; Hess , B. ; Lindahl , E. GROMACS: High Performance Molecular Simulations Through Multi-Level Parallelism from Laptops to Supercomputers . SoftwareX 2015, 1–2, 19–25 . (51). ↵ Abraham , M. J. ; van der Spoel , D. ; Lindahl , E. ; Hess , B. ; the GROMACS Development Team GROMACS User Manual. GROMACS Development Team , 2018 ; Version 2018. (52). ↵ Hornak , V. ; Abel , R. ; Okur , A. ; Strockbine , B. ; Roitberg , A. ; Simmerling , C. Comparison of Multiple Amber Force Fields and Development of Improved Protein Backbone Parameters . Proteins: Structure, Function, and Bioinformatics 2006 , 65 , 712 – 725 . OpenUrl (53). ↵ Petrova , S. S. ; Solov’ev , A. D. The Origin of the Method of Steepest Descent . Hist. Math . 1997 , 24 , 361 . OpenUrl (54). ↵ Lemak , A. S. ; Balabaev , N. K. On the Berendsen Thermostat . Mol. Simul . 1994 , 13 , 177 – 187 . OpenUrl CrossRef (55). ↵ Saito , H. ; Nagao , H. ; Nishikawa , K. ; Kinugawa , K. Molecular Collective Dynamics in Solid Para-Hydrogen and Ortho-Deuterium: The Parrinello–Rahman-Type Path Integral Centroid Molecular Dynamics Approach . J. Chem. Phys . 2003 , 119 , 953 – 963 . OpenUrl CrossRef (56). ↵ Petersen , H. G. Accuracy and Efficiency of the Particle Mesh Ewald Method . J. Chem. Phys . 1995 , 103 , 3668 – 3679 . OpenUrl CrossRef Web of Science (57). ↵ Hess , B. ; Bekker , H. ; Berendsen , H. J. C. ; Fraaije , J. G. E. M. LINCS: A Linear Constraint Solver for Molecular Simulations . J. Comput. Chem . 1997 , 18 , 1463 – 1472 . OpenUrl CrossRef PubMed Web of Science (58). ↵ Bjelkmar , P. ; Larsson , P. ; Cuendet , M. A. ; Hess , B. ; Lindahl , E. Implementation of the CHARMM Force Field in GROMACS: Analysis of Protein Stability Effects from Correction Maps, Virtual Interaction Sites, and Water Models . Journal of Chemical Theory and Computation 2010 , 6 , 459 – 466 . OpenUrl (59). ↵ Kovács , H. ; Mark , A. E. ; van Gunsteren , W. F. Solvent Structure at a Hydrophobic Protein Surface . Proteins: Struct., Funct., Bioinf . 1997 , 27 , 395 – 404 . OpenUrl (60). ↵ Chakrabarti , P. ; Pal , D. Main-Chain Conformational Features at Different Conformations of the Side-Chains in Proteins . Protein Eng . 1998 , 11 , 631 – 647 . OpenUrl CrossRef PubMed Web of Science (61). ↵ Chatterjee , P. ; Bagchi , S. ; Sengupta , N. The Non-uniform Early Structural Response of Globular Proteins to Cold Denaturing Conditions: A Case Study with Yfh1 . J. Chem. Phys . 2014 , 141 , 205103 . OpenUrl CrossRef PubMed (62). ↵ Hollingsworth , S. A. ; Karplus , P. A. A Fresh Look at the Ramachandran Plot and the Occurrence of Standard Structures in Proteins . BioMolecular Concepts 2010 , 1 , 271 – 283 . OpenUrl PubMed (63). ↵ Nelson , D. L. ; Lehninger , A. L. ; Cox , M. M. Lehninger Principles of Biochemistry , 5th ed.; W. H. Freeman and Company : New York , 2008 . (64). ↵ Jurrus , E. et al. Improvements to the APBS Biomolecular Solvation Software Suite . Protein Science 2018 , 27 , 112 – 128 . OpenUrl CrossRef PubMed (65). ↵ Kleywegt , G. J. ; Jones , T. A. Phi/Psi-chology: Ramachandran Revisited . Structure 1996 , 4 , 1395 – 1400 . OpenUrl CrossRef PubMed (66). ↵ Samanta , S. ; Chakrabarti , J. ; Bhattacharyya , D. Changes in Thermodynamic Properties of DNA Base Pairs in Protein–DNA Recognition . J. Biomol. Struct. Dyn . 2010 , 27 , 429 – 442 . OpenUrl PubMed (67). ↵ Mondal , M. ; Halder , S. ; Chakrabarti , J. ; Bhattacharyya , D. Hybrid Simulation Approach Incorporating Microscopic Interaction Along with Rigid Body Degrees of Freedom for Stacking Between Base Pairs . Biopolymers 2016 , 105 , 212 – 226 . OpenUrl PubMed View the discussion thread. Back to top Previous Next Posted May 05, 2026. Download PDF Supplementary Material Email Thank you for your interest in spreading the word about bioRxiv. NOTE: Your email address is requested solely to identify you as the sender of this article. Your Email * Your Name * Send To * Enter multiple addresses on separate lines or separate them with commas. You are going to email the following Simulation of Protein Structure using a Coarse-Grained Potential incorporating the Backbone Dihedral Interactions Message Subject (Your Name) has forwarded a page to you from bioRxiv Message Body (Your Name) thought you would like to see this page from the bioRxiv website. Your Personal Message CAPTCHA This question is for testing whether or not you are a human visitor and to prevent automated spam submissions. Share Simulation of Protein Structure using a Coarse-Grained Potential incorporating the Backbone Dihedral Interactions Kanika Kole , Abhik Ghosh Moulick , Jaydeb Chakrabarti bioRxiv 2025.08.20.671185; doi: https://doi.org/10.1101/2025.08.20.671185 Share This Article: Copy Citation Tools Simulation of Protein Structure using a Coarse-Grained Potential incorporating the Backbone Dihedral Interactions Kanika Kole , Abhik Ghosh Moulick , Jaydeb Chakrabarti bioRxiv 2025.08.20.671185; doi: https://doi.org/10.1101/2025.08.20.671185 Citation Manager Formats BibTeX Bookends EasyBib EndNote (tagged) EndNote 8 (xml) Medlars Mendeley Papers RefWorks Tagged Ref Manager RIS Zotero Tweet Widget Facebook Like Google Plus One Subject Area Biophysics Subject Areas All Articles Animal Behavior and Cognition (7636) Biochemistry (17704) Bioengineering (13898) Bioinformatics (41967) Biophysics (21460) Cancer Biology (18599) Cell Biology (25525) Clinical Trials (138) Developmental Biology (13384) Ecology (19909) Epidemiology (2067) Evolutionary Biology (24326) Genetics (15613) Genomics (22512) Immunology (17740) Microbiology (40423) Molecular Biology (17191) Neuroscience (88645) Paleontology (667) Pathology (2835) Pharmacology and Toxicology (4825) Physiology (7646) Plant Biology (15158) Scientific Communication and Education (2046) Synthetic Biology (4302) Systems Biology (9825) Zoology (2271)