Contents

On this page

Molecular Dynamics

Significance of molecular dynamics

Molecular Dynamics (MD) simulation is a fundamental computational methodology that uses statistical mechanics to evaluate the physical properties, structural stability, and mechanism of action of engineered proteins and peptides. By integrating Newtonian equations of motion, MD tracks the velocities, positions, and interactions between molecules over time. We applied this to our project to understand the translation and secretion efficiency of our protectants by analysing the interactions between different receptors.

With HADDOCK providing a structural framework for how our protectants bind to specific molecules, Molecular Dynamics (MD) simulations offer a glimpse into whether these predicted bound poses remain stable over time. Instead of treating a single docking pose as a definitive representation of the binding interaction, MD subjects the complex to a simulated dynamic environment, allowing the protein and protectant to move and interact as they would under molecular conditions.

Key findings and design implications

Results Histidine tag

The June and July design placed the purification tag at the C-terminus, which put six histidines beginning at Asn23 and the loss of the free carboxylate's participation in the native bidentate salt bridge with Arg487. While this tagged peptide satisfied the distance-based criterion used by the project, a more rigorous test considering both distance and bond angle showed that the bidentate salt bridge was missing in every simulated frame in contrast to 100% of native frames. The difference in these findings indicates the tagged peptide fails to maintain the native salt bridge, despite satisfying the project's original distance criteria. An N-terminal tag or synthetic peptide without a purification tag would not disrupt the C-terminal free carboxylate.

The same tagged construct under the project's own acceptance gate and under an atom-resolved criterion Six measures of one pair of trajectories. Under the three distance-only measures above the dashed line the tagged construct is indistinguishable from wild type and has the lowest conserved-window RMSD in the whole set, 1.12 angstroms. Under the three atom-resolved measures below it, every one reads 0 per cent against 99.8 to 100 per cent for wild type. One pair of trajectories, two definitions of the same contact bars are the measured value · 1,500 frames each wild type + 6×His BEFORE THE CONSTRUCT EXISTED the gate as written: distances between atoms Asn23-Arg487 within 4 Å either oxygen counts 100.0% 99.8% Clamp, 3 of 4 contacts held either oxygen counts 100.0% 99.7% Conserved-window RMSD lower is steadier 1.67 Å 1.12 Å AFTER, ON NAMED ATOMS the same trajectory, measured on the chemistry Bidentate salt bridge both oxygens, angle-aware 100.0% 0.0% Two or more H-bonds to Arg487 named donor and acceptor 99.8% 0.0% Carboxylate in Arg487 shell first solvation shell 100.0% 0.0%
Results pH 5.5

The slightly acidic apoplast environment led us to consider pH 5.5 instead of neutral pH as the default condition. At this pH, with His14, His22, and His227 protonated to reflect this environment, the bidentate salt bridge between peptide Glu12 and receptor His227 was present in 99.9% of frames compared to 0.8% at pH 7.0. This interaction is also present in the neutral pH AtPEP1 benchmark, and the acidic environment may facilitate a similar interaction in BoPEP4. We therefore found this to be the more suitable condition for studying BoPEP4, but we continued to simulate at pH 7.0 as a control.

Truncation crossed with protonation: only the combination fails A two-by-two. Full length at pH 7 is the reference. Protonating three histidines on the full-length peptide recovers the Glu12 to His227 contact from 0.8 to 99.9 per cent of frames and leaves the bidentate salt bridge at 99.9 per cent. The same protonation applied to the 9 to 23 construct drops bidentate occupancy from 46.9 to 6.9 per cent, the only run in the set to fail its acceptance criteria, and its backbone RMSD improves while it fails. The same two changes, crossed bidentate salt bridge (left) and Glu12-His227 (right), % of frames pH 7 pH 5.5, three histidines protonated full length 1-23 truncated 9-23 1-23 · pH 7 100.0 bidentate 0.8 Glu12-His227 reference RMSD 1.67 Å 1-23 · pH 5.5 99.9 bidentate 99.9 Glu12-His227 best untagged run RMSD 1.34 Å 9-23 · pH 7 46.9 bidentate 3.0 Glu12-His227 passed, but worst RMSD 2.28 Å 9-23 · pH 5.5 6.9 bidentate 26.4 Glu12-His227 failed the gate RMSD 1.33 Å Neither change alone pulls the C-terminus off the receptor: the His22 self-contact runs at 17.4% of frames from truncation and 64.7% from protonation. Only the two together reach 96.9%, and that is the frame count on which the clamp collapses. The AtPEP1 benchmark already holds Glu12-His227 at neutral pH, 3.36 Å in 76.9% of frames, so protonation recovers a contact the benchmark has.
Results Truncation ladder

Simulations found a salt bridge that the published hotspot lists we consulted do not carry: peptide Lys7 to receptor Asp129, held in 98.7% of frames in the 7–23 construct. Losing that anchor, and not the shortening itself, is what the 9–23 construct’s instability tracks, which is why stability does not fall step by step with each successive cut. The direct test is a full-length K7A mutant.

The truncation ladder is not monotonic: the four-rung comparison dips at 9 to 23 Clamp integrity and bidentate salt-bridge occupancy for the full-length peptide and the three cuts. Full length, 7 to 23 and 15 to 23 all hold the clamp above 98 per cent. Only the intermediate 9 to 23 construct falls, to 77.1 per cent clamp integrity and 46.9 per cent bidentate occupancy. The bottom row marks where the Lys7 to Asp129 salt bridge is present, at 96.6 per cent of frames for the full-length peptide and 98.7 per cent for 7 to 23, and where the residue has been cut away. Clamp integrity and the Arg487 salt bridge, by rung share of 1,500 frames · the ladder is not monotonic 0 50 100 gate 60% 100.0 100.0 1-23 full length 96.5 contacts RMSD 1.67 Å 98.5 87.7 7-23 Lys7 terminal 94.6 contacts RMSD 1.38 Å 77.1 46.9 9-23 ordered construct 68.5 contacts RMSD 2.28 Å 99.9 97.2 15-23 conserved core 89.6 contacts RMSD 1.81 Å clamp integrity, 3 of 4 native contacts bidentate Arg487 salt bridge Lys7-Asp129 96.6% 98.7% absent absent anchor, 2.85–2.81 Å
Results Wild-type BoPEP4 and the K18R mutation rationale

The wild-type run sent us back to the evidence for the planned K18R substitution. The Lys18 contact with Asp348 turned out to hold in only half the frames, across a wide swing in distance, and Phe371 spends most of the trajectory nearest Gly20 instead. The static docking model had overstated what Lys18 already does.

Simulation setup and parameters

All nine systems were solvated in explicit OPC water and simulated using the AMBER ff19SB force field with the PEPR1 ectodomain (residues 29-738) from PDB 5GR8. The systems have ~127,600 atoms each and underwent mass repartitioning of the hydrogen atoms to enable a 4 fs time step. Each trajectory ran 30 ns of production dynamics from a single random seed. We saved 1,500 frames from each one for analysis. This adds up to a total of 270 ns of simulation time or 25.9 GPU-hours on a single NVIDIA A100 GPU. The nine systems varied in four ways: the identity of the ligand, the length of the peptide, whether the peptide had a free C-end or not, and the protonation states of three histidines. All other simulation parameters were the same. The same analysis codebase was used for all nine trajectories, with identical definitions and calculations for all reported metrics, which allows for a common footing in analysis of structural stability and intermolecular interactions between the nine systems.

The nine systems, run by run
SystemLigandResiduesC-endHis+Clamp integrityBidentate Arg487RMSD Å
1–23reference BoPEP4 23 free 0 100.0 100.0 1.67
AtPEP1 7–23benchmark AtPEP1 17 free 0 99.9 99.9 1.17
7–23truncation BoPEP4 17 free 0 98.5 87.7 1.38
9–23truncation BoPEP4 15 free 0 77.1 46.9 2.28
15–23truncation BoPEP4 9 free 0 99.9 97.2 1.81
1–23 + 6×Histag BoPEP4 29 blocked 0 99.7 2.4 1.12
1–23pH 5.5 BoPEP4 23 free 3 100.0 99.9 1.34
9–23pH 5.5 BoPEP4 15 free 3 19.3 6.9 1.33
Apo receptorcontrol none no ligand · 0 · · ·

Wild type

The wild-type run served as a comparison and analysis benchmark for the other mutations. We found that the eight N-terminal residues had an average fluctuation above 2 Å. This was consistent with the alanine scan results, which found that mutations to the first six residues had the least penalty in the peptide-receptor docking, indicating a lower contribution to peptide-receptor interaction. All residues past Pro10 had over 75% contact with the receptor, except for Pro19. All residues past Arg11 had under 1 Å fluctuation, except for Lys18. Our conserved window 15-23, defined by aligning sequences from the PEP receptor family, had near 100% contact with mostly sub-Ångström movements. This evidence supports that the N- and C-termini have different behaviours and roles as the N-terminus has greater conformational freedom, while the conserved C-terminus region remains stabilised in contact with the receptor. The conserved residues also form the primary region of the peptide that interacts with the binding pocket, highlighting their importance.

The run also supported our decision to use truncations instead of point mutations. Arg9 has a fluctuation of 1.32 Å and a mean closest approach of 2.79 Å and 100% contact, which explains the sharp jump in the HADDOCK score when removing Arg9. Although Ile2 and Ile4 have contact with the receptor, they have a fluctuation of over 1.5 Å and switch between different locations on the PEPR1 instead of maintaining contact. Adding these residues also necessitates preserving the surrounding disordered residues, which have only about 25% contact with the receptor. The improved binding from adding these residues was not worth the extra sequence length, so we did not pursue the Δ1-2 and Δ1-3 designs and instead opted for the truncation designs with residues 7 and 9.

Per-residue contact occupancy and fluctuation for the wild-type reference run Twenty-three bars, one per residue. The opening GILIGS segment is low-occupancy and the fluctuation line sits above two angstroms across it. From Ser15 to Asn23 occupancy is near 100 per cent and fluctuation is sub-angstrom. The two ends of the peptide behave differently, and the conserved end is the one in contact. Contact occupancy (bars) and fluctuation (line), full-length BoPEP4 occupancy = share of 1,500 frames within 4 Å of the receptor conserved window 15-23 GILIGS, the cut 0 50 100 G 1 I L I G S K 7 K R 9 P R E P H S 15 S G K 18 P G G H N 23 2 Å fluctuation fluctuation, scale 0 to 4 Å Every residue past Pro10 holds over 75% contact except Pro19, and every residue past Arg11 fluctuates under 1 Å except Lys18. The two ends of this peptide do different jobs, and the end that can be cut is the end that is already loose.
Reassessing Arg18 contact network

The MD results also prompted us to reassess the rationale for the K18R substitution. The literature identifies two receptor contacts involving Arg18 in AtPEP1: one with Asp348 and another with Phe371. BoPEP4 contains lysine at the corresponding position, and our initial docking poses suggested that Lys18 could reproduce both interactions at approximately crystal-like distances. On this basis, we initially proposed K18R as a way to restore an interaction that was already largely established.

The MD trajectory, however, showed that these contacts were less stable than the docking poses suggested. Asp348 is the primary receptor partner of Lys18, with the next-closest peptide residue remaining 7.87 Å away, yet the Lys18–Asp348 contact is present in only 53.7% of frames, with distances ranging from 2.55 to 9.37 Å. By comparison, the Asn23–Arg487 interaction maintains 100% occupancy within a 0.54 Å distance window.

The interaction with Phe371 further illustrates the difference between the static model and the simulated behaviour. During the trajectory, Phe371 most frequently contacts Gly20, which remains within 3.69 Å for 84.4% of frames, while Lys18 averages 5.12 Å from Phe371 and approaches within 4 Å in only 18.1%. The initial docking poses overstated the extent to which Lys18 reproduces the interactions attributed to Arg18 in AtPEP1.

This does not rule out K18R as a viable design. Lysine has a shorter side chain than arginine, whose guanidinium group can support the proposed bidentate interaction. The distinction is that the substitution should not be described as simply restoring a contact that is already largely formed in BoPEP4. Introducing arginine is instead intended to provide the side-chain geometry needed to establish the proposed interaction more effectively. Whether K18R achieves this geometry remains untested and requires a dedicated simulation.

The stability of the Lys18–Asp348 contact also varies with peptide length. In the nine-residue construct, its occupancy increases from 53.7% to 90.1%, suggesting that the surrounding sequence influences how consistently the interaction forms. This difference highlights the importance of evaluating receptor contacts in the context of the complete peptide construct, as changes elsewhere in the sequence can affect the behaviour of individual residues.

Truncation ladder

An unidentified N-terminal anchor

The 9-23 construct was added to the truncation series to investigate the effects of N-terminal truncation on receptor engagement, and was the least stable of all the untagged constructs: 77.1% clamp integrity, the bidentate salt bridge in 46.9% of frames, the Asn23 to Arg487 distance inside the cutoff in 81.1%, and the largest backbone RMSD of the BoPEP4 trajectories at 2.28 Å. It formed 68.5 contacts in the conserved window against 96.5 in the wild type. It passed the acceptance gate, but every contact-based metric sat below the other untagged constructs run at pH 7. The salt bridge dissociated transiently instead of breaking once and staying broken, 94 times or 22.9% of the trajectory, with no break lasting longer than 0.68 ns, and Asp348 and Phe371 shifted their contacts between peptide residues. The interface was engaged, but much more dynamically than in the wild type.

We initially hypothesised that the unusual instability was caused by the novel N-terminus created by the deletion of the first 8 residues. In the 9-23 construct, the positively-charged alpha-ammonium group of Arg9 becomes the new N-terminal residue, which could be forming a new interaction with the receptor. To investigate this possibility, we simulated the 7-23 construct, which moves the free N-terminus to Lys7, which is an internal residue. The interface became much more stable, the clamp integrity jumped to 98.5% from 77.1%, the occupancy of the Asn23-Arg487 interaction rose to 99.9% from 81.1%, and the occupancy of the bidentate salt bridge jumped to 87.7% from 46.9%. The shared 9-23 sequence had its contacts increasing to 151.3 from 128.7.

These changes supported the importance of the residues removed in the 9-23 construct, but did not prove that it was due to the charge on the N-terminus. In the 7-23 construct, the alpha-amino group of Lys7 was 7.11 Å away from its closest receptor carboxylate and within 4 Å of it in 0.0% of frames. In the 9-23 construct, the N-terminal Arg9 was within 4 Å of a receptor carboxylate in only 18.5% of frames. The full-length BoPEP4 and the AtPEP1 benchmark had little evidence of their N-termini making persistent interactions.

The truncation series had revealed an unexpected N-terminal anchor that was not included in the published hotspot list. The side chain of Lys7 formed a salt bridge with Asp129 at 2.81 Å in 98.7% of frames, even when it was the N-terminal residue in the 7-23 construct. Arg9 contributed to the interaction network in 100% of frames. These identified two N-terminal salt bridges in the BoPEP4 that were not in the published hotspot list and offered an alternative explanation to the instability of the 9-23 construct that did not rely on the introduction of the interaction with the terminal alpha-ammonium group.

The instability of the 9-23 construct is more consistent with the removal of the Lys7-Asp129 anchor while residues 9-14 were retained within the binding interface, which is why it is not consistently less stable with decreasing peptide length: full-length BoPEP4, 7-23 and 15-23 all hold clamp integrity above 98%, and only the intermediate 9-23 drops to 77.1%. What matters is which interactions were removed, not how many residues remained. The direct test is a full-length K7A mutant, which would remove the proposed side-chain anchor without changing the peptide termini.

7–23 as a benchmark for AtPEP1

The 7-23 construct offers a more direct comparison to the AtPEP1 crystal structure as they both span the same 7-23 sequence and share the same free N-terminus. 15 out of the 17 positions had very similar contact occupancies between the two peptides. The main sequence difference within the published receptor-contact hotspots was that AtPEP1 had an Arg18 while BoPEP4 had a Lys18, which is the structural basis for the proposed K18R substitution.

There was a difference in the anchoring of the N-termini between the two peptides. AtPEP1 used Gln8 to make a contact with Tyr153 at 3.41 Å in 95.8% of frames while Lys7 contributed in only 19.0% of frames. Conversely, BoPEP4 used Lys7 to make a contact with Asp129 at 2.81 Å in 98.7% of frames while Lys8 contributed in only 4.1% of frames. One residue was primarily making the contact at this end of the interface for both peptides, but the anchoring residue and receptor partner were different by both sequence and position. They shared a free N-terminal residue of Lys7 and this comparison is less affected by the difference in the terminal configuration and more directly compares the different interaction networks of the two peptides. This is also seen in the occupancy of the Lys7 contact, which reached 99.9% in the BoPEP4 and only 23.5% in the AtPEP1.

Otherwise, 7-23 closely recapitulated the binding architecture of full-length BoPEP4, with two exceptions: the Asn465 hydrogen bond held in 62.9% of frames against 99.8% in wild type, and the bidentate salt bridge in 87.7% against 100%. It carries all the main features of the interface at reduced persistence, so it benchmarks the terminal configuration and does not stand in for the full-length interaction network.

15–23 as a negative control

The 15-23 construct was designed to be a negative control as the sequence before it was proposed to be important for stable receptor engagement. Unexpectedly, the peptide was highly organised throughout the simulation. It had a clamp integrity of 99.9%, occupancy of the bidentate salt bridge of 97.2%, and a backbone RMSD of 1.81 Å. Its nine residues formed 89.6 contacts compared to 96.5 contacts made by the same residues in wild type. Despite the deletion of 14 residues, the construct had ~93% of the contact network of its corresponding sequence in the full-length peptide. However, the dissociation of peptides from a leucine-rich repeat (LRR) ectodomain could have microsecond to millisecond timescales, much longer than the 30 ns trajectories simulated here. The lack of dissociation in the simulation does not prove that the peptide would not eventually dissociate or that its equilibrium affinity is unchanged. Still, it shows that the conserved C-terminal region can maintain an organised binding register for the duration of the simulated trajectory in the absence of the sequence before it.

Implications of the truncation series

The two constructs that bracket 9-23 make the point. 15-23 lacks the intervening 9-14 sequence altogether and still maintains an organised interaction through its conserved C-terminal contacts, while 9-23 keeps that sequence and loses the anchor holding it.

The truncation series indicates that the conserved C-terminal region forms the structural basis for receptor engagement, while the sequence before it contributes extra interactions that affect the stability and dynamics of the interface.

pH 5.5

Because conventional molecular dynamics does not dynamically simulate pH, our pH 5.5 condition represents a fixed protonation state intended to approximate the apoplast, not a constant-pH simulation. Histidine protonation was assigned using model pKa values, with His14 and His22 on BoPEP4 and His227 on PEPR1 carrying positive charges, while all other histidines remained neutral and all Asp and Glu residues remained deprotonated throughout the 30 ns trajectory. No other simulation parameter was changed relative to the wild-type pH 7 run, making the three protonation changes the sole variable in this comparison.

Under this condition, the interaction between peptide Glu12 and receptor His227 was substantially strengthened, with the pair separated by 2.75 ± 0.11 Å in 99.9% of frames compared with 5.33 Å and only 0.8% occupancy at pH 7. This interaction is already present in the AtPEP1 benchmark at neutral pH, where the corresponding pair occurs at 3.36 Å with 76.9% occupancy, indicating that protonation allows BoPEP4 to recover a contact present in its structural benchmark instead of introducing an entirely new interaction. This improvement did not compromise the primary C-terminal anchor: clamp integrity remained at 100.0%, the bidentate salt bridge at 99.9%, and the conserved window formed 95.7 contacts with a backbone RMSD of 1.34 Å, the lowest value observed among the untagged constructs that satisfied their acceptance criteria.

This result therefore changes the interpretation of the original sequence difference at Glu12-His227. The lower contact occupancy observed at pH 7 had initially appeared to reflect a sequence-level deficiency in BoPEP4, but the pH-controlled comparison indicates that it was instead a consequence of the protonation state used in the simulation. Because the apoplast is the environment in which the engineered protectant is expected to function, pH 5.5 is therefore a more physiologically relevant default for full-length BoPEP4, with pH 7 retained as the control. At the same time, the trajectory does not support the proposed His14 pH-switch mechanism: His14 remained approximately 6.1 Å from the nearest receptor carboxylate and formed a sub-4 Å contact in 0.0% of frames at either condition. The pH-responsive interaction identified here is instead the Glu12-His227 pair.

Evaluating protonation in the 9–23 construct

Protonation improved the interaction network of full-length BoPEP4 at pH 5.5, raising the question of whether the same effect could improve the stability of the 9–23 construct, which was the least stable untagged peptide at pH 7. To investigate this, we repeated the 9–23 simulation using the same starting coordinates and protonating the same three histidines as in the full-length pH 5.5 simulation, while keeping all other simulation parameters unchanged.

Contrary to the full-length result, protonation further destabilised the 9-23 construct. Asn23-Arg487 remained within 4 Å in only 23.6% of frames, well below the 60% acceptance threshold, while clamp integrity fell to 19.3% and bidentate salt-bridge occupancy to 6.9%, making this the only BoPEP4 trajectory that failed its pre-registered acceptance criteria. The underlying mechanism was a new intramolecular interaction between protonated His22 and the C-terminal carboxylate: this contact occurred at 2.84 Å in 96.9% of frames in 9-23 at pH 5.5, compared with 17.4% at pH 7. For comparison, the corresponding interaction occurred in only 1.1% of full-length frames at pH 7 and increased to 64.7% at pH 5.5 while the full-length clamp remained intact.

Neither condition alone was sufficient to produce the strong self-contact observed in the failed trajectory: truncation increased occupancy from 1.1% to 17.4% at neutral pH, while protonation increased full-length occupancy to 64.7%; only their combination produced the 96.9% occupancy that coincided with collapse of the C-terminal anchor. At the same time, the stabilising Arg9-Asp179 interaction remained highly occupied in full-length trajectories but decreased to 75.0% in 9-23 at pH 5.5, leaving the shortened peptide sufficiently unconstrained for its C terminus to fold back toward His22. Protonation should therefore not be treated as a universal stabilising condition across constructs; instead, its effect depends on the interaction network remaining in each peptide.

Two additional observations support this interpretation. First, the newly exposed Arg9 α-amino group remained 6.90 Å from its nearest receptor carboxylate and formed a sub-4 Å contact in only 5.8% of frames, providing further evidence against the proposed terminal salt bridge. Second, although protonation restored the Glu12-His227 interaction to 99.9% occupancy in full length, the same pair reached only 26.4% in 9-23, indicating that protonation alone cannot stabilise an interaction when the upstream peptide sequence required to position it has been removed.

Apo receptor

To determine whether peptide binding itself produced measurable conformational changes in PEPR1, we also simulated the receptor without peptide using chain A from the wild-type complex and otherwise identical conditions. Contrary to the expectation that the unbound receptor would become substantially more flexible or that its binding groove would partially collapse, the apo receptor remained structurally comparable to the peptide-bound trajectories.

The apo receptor had the lowest mean Cα RMSD of all nine simulations at 1.57 Å, compared with 1.58-2.16 Å across the peptide-bound runs, while its groove span was 35.91 Å compared with 35.65-36.42 Å for the bound structures. Its per-residue flexibility profile also correlated strongly with the wild-type profile at r = 0.94-0.97. The 54 interface residues shifted by only +0.125 Å relative to wild type, compared with +0.156 Å for the remaining 656 residues, with the largest differences occurring at residues 41-55, which constitute the N-terminal cap of the ectodomain, and with residue 738, the final residue, being the most rigid in the structure. Both observations reflect chain-end effects distant from the binding groove. Within the 30 ns simulation window, we found no evidence for a substantial peptide-induced conformational change in the receptor.

This provides both a structural and methodological control. The pre-formed groove supports representing the binding surface using a mean receptor structure instead of assuming that ligand binding creates the pocket. More importantly, the apo trajectory establishes a noise floor against which receptor-side differences can be interpreted: the seven peptide-bound runs differed from one another by 0.59 Å in receptor Cα RMSD, while the apo trajectory was only 0.01 Å below the steadiest bound trajectory. Because changing the peptide produced more receptor-side variation than removing it entirely, we do not interpret small receptor RMSD or flexibility differences as evidence that BoPEP4 stabilises PEPR1. Instead, our comparisons remain centered on peptide-side contact geometry.

Receptor alpha-carbon RMSD across all nine runs, with the apo run as the floor Nine rows. The receptor simulated with no peptide at all has the lowest RMSD of the set at 1.57 angstroms, one hundredth of an angstrom below the steadiest bound run. The bound runs spread over 0.59 angstroms among themselves. Changing the peptide moves the receptor more than removing it, so no receptor-side difference in this set is read as peptide-induced stabilisation. Receptor Cα RMSD, all nine runs the receptor on its own is the steadiest trajectory in the set apo, no peptide 1.57 1-23 + 6xHis 1.59 1-23 1.58 7-23 1.63 AtPEP1 7-23 1.70 9-23 1.72 1-23 pH 5.5 1.87 9-23 pH 5.5 2.14 15-23 2.16 apo floor 1.57 Å 0.59 Å separates the bound runs from each other; 0.01 Å separates the apo run from the steadiest of them. Changing the peptide moves the receptor more than removing it, so no receptor-side difference here is read as stabilisation.

Overall

The combined trajectories demonstrate why RMSD and aggregate structural measures are insufficient for ranking these constructs. Across the eight peptide-containing runs, backbone RMSD over the conserved window correlated only weakly with the bidentate Arg487 salt bridge, with Pearson and Spearman coefficients of +0.13 and +0.29, respectively, while correlations with clamp integrity were −0.03 and −0.07.

The individual trajectories illustrate this more clearly: the tagged construct and 9-23 at pH 5.5 had the two lowest RMSDs, at 1.12 Å and 1.33 Å, yet their bidentate salt bridges were occupied in only 2.4% and 6.9% of frames, whereas 15-23 had a higher RMSD of 1.81 Å but retained the salt bridge in 97.2% of frames. Contact count performed relatively better, with a Spearman correlation of +0.71 across the seven untagged runs, but this comparison is based on only seven trajectories and is strongly influenced by individual constructs. Together with the earlier docking-score analysis, these results support using named atom-pair criteria and not aggregate scores to determine whether a construct retains the intended binding chemistry.

Conserved-window backbone RMSD against bidentate salt-bridge occupancy, eight runs A scatter of the eight peptide-bearing runs. The horizontal axis is backbone RMSD over the conserved window, the vertical axis is the share of frames holding the bidentate Arg487 salt bridge. The points do not line up: the tagged construct at 1.12 angstroms and the failed 9 to 23 run at pH 5.5 at 1.33 angstroms are the two lowest RMSDs in the set and hold the salt bridge in 2.4 and 6.9 per cent of frames, while the 15 to 23 control at 1.81 angstroms holds it in 97.2 per cent. Backbone RMSD does not predict whether the salt bridge survives eight peptide-bearing runs · Pearson +0.13, Spearman +0.29 vertical: bidentate Arg487 salt bridge, % of 1,500 frames 0 25 50 75 100 1.2 1.6 2.0 2.4 conserved-window backbone RMSD (Å) · steadier to the left AtPEP1 7-23 crystal benchmark 1-23 reference 1-23 pH 5.5 best untagged run 15-23 control that held 7-23 9-23 ordered construct 9-23 pH 5.5 failed the gate 1-23 + 6xHis tag The two steadiest runs in the set, by the number a construct would normally be ranked on, are the tagged construct and the only run that failed its gate. Contact count does better, Spearman +0.71 over the seven untagged runs, but on seven points. Named atom pairs do better than either.

Limitations

The 30 ns timescale limits what can be inferred about binding strength. Protein-peptide dissociation happens on timescales orders of magnitude longer, so a peptide that stays in the groove for a whole trajectory has shown structural persistence over that window and not strong binding. Settling binding strength needs longer runs, independent sampling, and free-energy calculations or an experimental measurement.

The absence of replicate simulations limits how far the construct comparisons can be trusted. Each system ran as a single trajectory, so a single set of initial velocities and one relaxation pathway, and systems like these can relax along different pathways even from similar starting structures under identical parameters. Differences in contact occupancy, backbone fluctuation or structural integrity could therefore be variance and not a consequence of peptide design. Multiple independent seeds per system would settle whether the observed behaviours generalise and whether the differences between constructs are robust.

The acceptance criteria need care when the tagged and untagged constructs are compared. They were set before the C-terminal histidine-tagged construct was simulated, around interactions present in the untagged peptide, chiefly the one between the free C-terminal carboxylate of Asn23 and Arg487 on PEPR1. Adding six histidines after Asn23 changes that chemistry: Asn23 is no longer the free C-terminus, so the original interaction cannot be assessed in the same way. The trajectory is still usable, but the criteria should be revised to separate the interactions that survive the tag from those that depend on the native free C-terminus.

The acidic-pH simulations are a specific model, not a comprehensive account of pH-dependent behaviour. His14 and His22 on BoPEP4 and His227 on PEPR1 were given a positive charge from their estimated pKa, every other histidine was left neutral, and all of it was held fixed for the trajectory, so no residue could gain or lose a proton as its local environment changed. The Glu12-His227 interaction is therefore a pH-dependent interaction under these conditions and not a demonstration of how it behaves across the full range found in the apoplast. The constraint bites hardest on the tagged construct, whose histidines were modelled as neutral even under the acidic condition, so this run does not show how the tag itself behaves there. A constant-pH scheme, or better-chosen tag protonation states, would resolve it.

Application of results

Wetlab

The protein design workflow chart connects the computational findings to the experiments that follow because simulations can help identify which design changes are worth testing, but their predicted stability must ultimately be evaluated through biological experiments. For the histidine-tagged construct, the results informed the decision to move the purification tag to the N-terminus before expression, preserving the native C-terminal chemistry implicated in receptor binding. Similarly, the relatively consistent contact network across the longer and shorter constructs raised the possibility that truncation could reduce the peptide's production burden without eliminating its key interactions. This remains a manufacturability hypothesis, however, and will be tested through secretion assays comparing the truncated constructs with wild-type BoPEP4. These experiments will determine whether the computationally identified designs can be produced and secreted effectively in practice.

References

  1. Tian, C. et al. ff19SB: amino-acid-specific protein backbone parameters trained against quantum mechanics energy surfaces in solution. Journal of Chemical Theory and Computation 16, 528 to 552 (2020).
  2. Izadi, S., Anandakrishnan, R. & Onufriev, A. V. Building water models: a different approach. Journal of Physical Chemistry Letters 5, 3863 to 3871 (2014). OPC.
  3. Tang, J., Han, Z., Sun, Y., Zhang, H., Gong, X. & Chai, J. Structural basis for recognition of an endogenous peptide by the plant receptor kinase PEPR1. Cell Research 25, 110 to 120 (2015). PDB 5GR8.