Computational Peptide Design Cotargeting Glucagon and Glucagon-like Peptide-1 Receptors

Peptides are sustainable alternatives to conventional therapeutics for G protein-coupled receptor (GPCR) linked disorders, promising biocompatible and tailorable next-generation therapeutics for metabolic disorders including type-2 diabetes, as agonists of the glucagon receptor (GCGR) and the glucagon-like peptide-1 receptor (GLP-1R). However, single agonist peptides activating GLP-1R to stimulate insulin secretion also suppress obesity-linked glucagon release. Hence, bioactive peptides cotargeting GCGR and GLP-1R may remediate the blood glucose and fatty acid metabolism imbalance, tackling both diabetes and obesity to supersede current monoagonist therapy. Here, we design and model optimized peptide sequences starting from peptide sequences derived from earlier phage-displayed library screening, identifying those with predicted molecular binding profiles for dual agonism of GCGR and GLP-1R. We derive design rules from extensive molecular dynamics simulations based on peptide–receptor binding. Our newly designed coagonist peptide exhibits improved predicted coupled binding affinity for GCGR and GLP-1R relative to endogenous ligands and could in the future be tested experimentally, which may provide superior glycemic and weight loss control.


S1. Convergence of MD simulations
To evaluate the convergence of MD simulations, we monitored the secondary structure elements of receptors throughout the 100 ns MD runs. We computed the timelines of cumulative average secondary structure features. Cumulative averages calculate the average property at every timestamp of simulation, such that each sliding averaged property over the time interval is computed. Figs. S4 and S5 show the preservation of protein secondary structures in the GCG and GLP-1 receptor throughout the simulations and provide support for the high stability of the receptor structure in the complex system between co-agonist and GRs. Fig. S6 depicts the low Root Mean Square Fluctuation (RMSF) of GCGR and GLP-1R, respectively during the simulations with enumerated agonist peptides. RMSF, i.e. standard deviation of atomic position in the 100 ns trajectory after fitting to the reference "average" structure of peptide/receptor, is calculated to check the fluctuations of C-alpha atoms of receptor in individual co-agonist/GR complex in the simulations. The highly fluctuating atoms in residues in all the systems (including the apo-GRs) are in the intra-and extracellular loop regions while other regions displayed low fluctuation over time (Fig. S7). Since the GRs ECD consists of many mobile loops, it exhibits a larger deviation than the transmembrane domain (TMD) for both receptors. The root-mean-square deviation (RMSD) analysis allows the quantification of the degree of conformation changes that occurred during the co-agonist/GR simulations. To examine the conformational changes in GRs, we used the post-processed MD trajectories of the receptors to generate the RMSD plots (Fig. S8). It was noted that the GCGR backbone RMSDs exhibited deviation in the range of 0.5-2.5 Å concerning their equilibrated structure, mostly this deviation was caused due to the structural changes in flexible termini and loop regions. The intrinsic dynamics stability of the designed PDL dual co-agonist against the endogenous peptide ligand (GCG and GLP-1) during the simulation has been summarized in Table S10, which shows the peptide helix is very stable throughout the simulation. The peptide remains stable and intact during the simulation, as formulation development of parenteral peptide therapeutics frequently encounters aggregation challenges. Native contacts determine protein folding and stability mechanisms in atomistic simulations. It is evident from the fraction of native contacts analysis, that both the receptor in the simulations with co-agonistic peptides, reached a plateau after 50 ns (except for P11/GLP-1R complex which reflects the stability after reaching ~80 ns) from the fraction of native contacts (Fig. S9).
Fraction of native contacts, Q(X) 1 presents the property of retainment of contacts in the native folds of receptors as a function of simulation time, and was calculated using the following equation: Where, N represents the set of all pairs of heavy atoms (i, j) that are in contact if their distance is less than 6 Å and are separated by at least 3 residues. The distance between heavy atoms i and j in the conformation (x) sampled at time t is represented by rij(x). The smoothing parameter, β, is taken to be 5 Å -1 , and λ is a factor that describes fluctuations when the contact is formed, with a value of 1.8, which are standard values 1 .

S2. Calculation of binding energy
Exploring co-agonist and GRs interaction with MD simulations and binding free energy calculations: Quantitative knowledge of class B1 GPCRs (glucagon and glucagon-like peptide-1 receptor)-peptide ligand binding affinities is essential in understanding molecular recognition; hence, efficient, and The free energy term Gx, where x corresponds to GPCR/co-agonist or GPCR or c-agonist and is given by, where Ebonded comprises bond-stretch, angle-bend, torsion, and improper-dihedral energies, and EvdW and Eele are the van der Waals and electrostatic nonbonded interaction energies, respectively.
Together, the sum of these terms makes up the vacuum molecular mechanics (MM) energy terms as shown below, while Gpolar and Gapolar in Equation S3 constitute the solvation-free energies calculated using a continuum (implicit) solvation model, representing the free energy change due to transferring a solute in a vacuum to the solution. Gpolar is the electrostatic contribution to solvation and is obtained by solving the Poisson-Boltzmann (PB) equation 2 , while the Gapolar term is the nonpolar contribution and S6 is often approximated by a solvent accessible surface area (SASA) term. Gapolar is estimated from a linear relation to solvent-accessible surface area (SASA) as: where γ is a coefficient set to the surface tension of the solvent and b is a fitting parameter. TS is the entropy contribution, where T is the absolute temperature and S is the configurational entropy. It should be noted that although the solvation energy change, ∆G solv is a free energy term, the ∆E MM does not consider the change in entropy due to binding, and hence is not a free energy term. Therefore, the total binding energy ∆G bind approximates the Gibbs free energy, with entropy contributions to the overall binding free energy difference for one peptide ligand vs. another, the binding specificity, ∆∆G binding considered minimal in the present case as both ligands both bind to the same site and are similar in size, typical of competitive native vs. inhibitor ligand binding to biological receptors.
Free energy calculation methods including MM/PBSA or MM/GBSA are suitable for predicting the binding free energies and evaluating the relative stabilities and binding strengths of different biomolecular complexes with a low computational cost 3,4 . We note that the slight overall unfavorability observed in peptide-receptor affinities for a few systems (see Table S2) stem from larger penalties imposed by polar solvation energies (solute-water electrostatics) on solute-solute electrostatic energies estimated in a vacuum. It is important to note that a positive overall G does not necessarily reflect an unfavourable binding of the peptide to the receptor 5 . The calculated binding free energy is only an estimate of the change in free energy associated with the process of ligand binding to the receptor and could be accurately evaluated from their relative binding energies (G) highlighting the specificity of peptide binding in one system over the other (see Fig. S10D).
Additionally, there may be other intrinsic factors such as entropic contributions and changes in the solvent environment that could affect the overall binding affinity 6 .

S3. Computed interaction and contact maps
Interaction energy profiles were calculated using GROMACS tools 7 . Interdomain contact probability maps and interaction maps were generated using the CONAN contact analysis tool 8 . These interaction maps identify H-bonds, salt bridges and hydrophobic networks between the co-agonist peptide ligand and GCG/GLP-1 receptors. These maps were produced with a cut-off inter-atomic (heavy atoms) distance of 5 Å as the minimum average distance between the heavy atoms of each residue. The intermolecular interaction maps were then generated using a truncation lifetime of 0.5. The average S7 numbers of hydrogen bonds (H-bonds) were counted using the auxiliary GROMACS 7 module gmx hbond using the criteria of donor-acceptor distance ≤ 3.5 Å and hydrogen-donor-acceptor angle ≤ 30°. The choice of 3.5 Å cutoff distance includes also weaker H-bond interactions between 3.0 and 3.5 Å, which play a crucial role in determining protein secondary structure conformations 9, 10 .

S4. Free energy maps
We computed 2D free energy maps of GR ligand-bound and unbound states using Root Mean Square Deviation (RMSD) of backbone atoms and radii of gyration (Rgyr) order parameters (Fig. S26), and ligand-receptor interaction energy and receptor backbone Root Mean Square Fluctuations (RMSF) (Fig. 27). The values for the Gibbs free energy (kJ/mol) using two order parameters were obtained using the equation: where Pi is the probability distribution for pairs of order parameters, and Pmax is its maximum, such that lnPi -lnPmax identifies the lowest free energy point at ΔG = 0.

S5. Peptide point mutation in the design of co-agonists
Five selected PDL-peptide sequences were modelled using the experimental GCG helical peptide structure as a reference template (see Methods) to construct the PDL-peptides : GR complexes (P11, P23, P28, P32, and P35 each in complex with GCGR and GLP-1R) and the "MD-directed design" (MDD) peptides in complex with GRs (see Table S5 for mutation points on the peptide templates) as a starting point for MD runs. Full-length GCG : GCGR and GLP-1 : GLP-1R complexes were also modelled to guide systematic mutagenesis and rationally improve co-agonist selectivity to GRs. Nterminal residues (His1 and Ser2 or Ala2, for GCGR and GLP-1R, respectively) of peptides are known to be critical for GR interaction 11 and potency 12 and mutations at these points of agonist peptides may lead to antagonistic effect 13 . Previous studies have reported that substitution at the 22 nd or 23 rd position in the C-terminus of the peptides leads to a lower affinity for GRs which lowers receptor activation [14][15][16] . Hence, we retained these N-terminus residues when designing the GR coagonists.
Affinity and residue specificity to GRs of the designed point mutations on the PDL-peptides was compared against the native peptide endogenous ligands (see Tables S6, S7) by accounting for the contributions at all residue positions towards overall net ΔGbind. Q20H substitution on P11 and P32 S8 favoured binding to GLP-1R as position 20 also increases the affinity of native GLP-1 to GLP-1R.
The effect is analogous to the positively charged side chain of Lys20 in GLP-1 facilitating receptor binding 17 .

S6. MD-directed Design (MDD) of a GR co-agonist
Based on the data obtained from PDL peptide/GR molecular level simulations along with endogenous ligands, we designed and tested three new GR co-agonist peptides (MDDGCGR, MDDGLP-1R, MDDGR).

S7.1 Insights into the structural basis of GR activation by peptide co-agonists
Both GRs, GCGR and GLP-1R send activation signals mainly via the G protein subtype (Gs) class of heterotrimeric G proteins [18][19][20] . While a recent study 20 reports that the GCGR activation may present some aberration in G protein coupling by also binding to a different G protein subtype (Gi), in both cases, the glucagon initially binds to the same extracellular site. Briefly, extracellular binding to class B1 GPCRs triggers a conformational change that initiates heterotrimeric G protein coupling on the intracellular receptor domain. This stimulates the GPCR to transform the inactive G protein into its active form, which further enables the intracellular physiological processes. Our computational models also suggest a large degree of GR conformational dynamics that is required for the co-agonist binding and receptor activation. We note that the designed PDL and MDD co-agonist occupies the native active site cleft within the GR's 7-TM helical bundle (see Figs. S27B, G).