Identification and engineering of highly functional potyviral proteases in cells using co-evolutionary models

Domain: chemistry
Published: 2026-02-26
Editorial outcome: minor_revision

Abstract
Efficiency and substrate specificity of proteases in the Potyviridae family have not been comprehensively profiled. Here we develop a model that learns co-evolutionary features to accurately predict and experimentally validate protease performance at single amino-acid resolution. We identify and engineer several proteases that perform better than the commercially available tobacco etch virus protease. To demonstrate the resolving power of our methods, we engineer protease crosstalk to selectively trigger a synthetic cell-death program in human cells.

Introduction
Proteases are commonly used tools in protein chemistry1,2, industrial applications3, and as composable parts of synthetic biological circuits4–7. Nuclear Inclusion a (NIa) proteases are cysteine proteases in the Potyviridae family that can recognize a unique cleavable sequence (cs) of seven amino acids8 (Fig. 1a,b). For example, the tobacco etch virus protease (TEVp) is known to cleave its TEVcs: ENLYFQ^S (where ^ is the cleaved peptide bond) and is commercially produced for use in biochemical purification of recombinant proteins. More than 3800 NIa proteases and variants are known, but their catalytic activities and substrate specificities have yet to be determined.

We hypothesized that co-evolutionary features of protease-substrate sequences can be used to quantitatively predict and validate the function of Potyviral proteases in human cells. The co-evolutionary pressure to maintain protease-substrate specificity influences which pairs possess the sequence composition to perform cleavage and therefore persist in nature. To investigate the specificity of protease-substrate interactions, we have compiled 3817 pairs of NIa protease and substrate (Supplementary Data 1) that provide a statistical setting to build the Protease Substrate Specificity Calculator (ProSSpeC). ProSSpeC leverages direct coupling analysis (DCA) to learn patterns of specificity by estimating a joint probability distribution of sequences in interacting families. DCA has been used to study protein structure 9, complex formation 10, conformational plasticity 11 and specificity between two protein families 12,13. It also served as inspiration for the Evoformer module of AlphaFold2 and the Pairformer module of AlphaFold3 14,15, which have also been utilized for protein ligand docking16. DCA can identify epistatic interactions, and we use the collection of DCA parameters learned from paired protease-substrate sequences to construct a Hamiltonian specificity (H ) score (Fig. 1c). Any change in sequence composition (of both protease and substrate) affects the strength of amino acid couplings and the resulting H score, providing a quantitative measure of specificity upon mutation. We reasoned that ProSSpeC can be used to predict and design protease sequences that can cleave specific substrate sequences.

Here we show how our quantitative modeling method, ProSSpeC, can resolve single amino-acid effects on protease-substrate interaction. We first establish that ProSSpeC can predict proteolytic cleavage efficiency against any natural potyviral substrate sequence. We then show that ProSSpeC can predict changes in cleavage efficiency of substrates undergoing point mutations. Last, we apply the predictive power of ProSSpeC to engineer selective protease-induced cell death at single amino-acid resolution. Together, these results show that ProSSpeC can accelerate the design of proteases for various biomedical and biotechnology applications.

Results
ProSSpeC is a coevolutionary model that leverages the covariation of NIa protease sequences along with their aligned substrate sequences to calculate a Hamiltonian specificity score (H ), where the more negative the H score is, the stronger the predicted specificity is. This specificity score is defined by the Potts Hamiltonian energy of canonical protease-substrate pairs (Equation 1), with the Potts Hamiltonian energy of scrambled enzyme-substrate pairs (Equation 2) removed, to provide a score of sequence attributes unique to interacting pairs and not shared due to family-likeness or general sequence similarity (Equation 3). For both canonical and scrambled protease-substrate pairs, the Potts Hamiltonian is defined by as the sum of all parameters in a Boltzmann distribution, with couplings (e ) and local fields (h) are inferred using direct coupling analysis (DCA), specifically the mean-field formulation (mfDCA) 9. A similar rationale for using the DCA Hamiltonian to characterize molecular interactions has been applied to study specificity in two-component systems12, recognition in protein-RNA interactions17, and to predict and engineer compatibility in hybrid transcription factors13,18. The ProSSpeC model applies this approach to quantify interactions of the protease and the substrate. When we were interested in only the 7-amino acid substrate motif; we masked out substrate residues outside of this window in H calculations to yield the masked H (Equation 4). The masked H was only utilized for predictions where the designed experiment did not include the extended substrate context. The ProSSpeC model allows us to quantitatively assess the specificity of NIa proteases and their substrate motif through sequence alone, allowing us to filter and rank both natural NIa protease-substrate pairs as well as mutated ones prior to experimental testing.

To rapidly screen for sequence-specific cleavage activity, we repurposed a fluorescent protein reconstitution assay in human cells that has been used to cleave and release protein cargoes into the nucleus19 (Fig. 1d). We initially assayed a set of 31 protease-substrate sequence pairs (Supplementary Fig. 1), which included a few proteases with favorable (i.e., most negative) ProSSpeC H scores (Fig. 2a). A subset of these proteases exhibited cleavage-induced fluorescence that are equal to or better than that of TEVp, and they also possessed similar or more favorable H scores than TEVp (Fig. 2b). This result supports our hypothesis that ProSSpeC can infer proteolytic cleavage activity from co-evolutionary sequence features of protease-substrate families.

To further determine if ProSSpeC can estimate protease specificity against any substrate sequence, we determined the specificity of fifteen proteases against fifteen substrates in human cells. We found that favorable H scores arise when a protease is paired with its putative native substrate sequence (Fig. 3a). For example, when pairing TEVp with each of the fifteen substrates, the TEVp-TEVcs pair has the lowest H score (H = -147; H-range of [-147, 1]) and exhibited the highest cleavage-induced fluorescence intensity than any other TEVp-substrate pairs. Importantly, the fluorescence intensities of all 225 protease-substrate pairs are correlated with H scores (Fig. 3b; see also Supplementary Fig. 2). Also, using these pairs and to quantify how the model captures crosstalk in off-diagonal elements, we generated a Receiver Operating Characteristic (ROC) curve with an Area Under the Curve (AUC) = 0.878 showing that the model exhibits strong discriminatory ability between true cleaving and non-cleaving pairs (H cutoff = -84) (Supplementary Fig. 3). These results demonstrate the predictive power of ProSSpeC and our ability to experimentally validate protease-substrate specificity in human cells.

Our experimental assay also revealed a subset of proteases that exhibited minimal or no cleavage-induced fluorescence (Supplementary Fig. 1). Two potential scenarios could explain the apparent lack of cleavage: (1) protease is autoinhibited, which has been shown previously for wild-type TEVp20; or (2) the putative substrate sequence does not promote cleavage18. We addressed the first scenario by mutating the residue position known to be involved in autoinhibition20 (Supplementary Table 1). To address the second scenario, we retested nine proteases that exhibited minimal or intermediate cleavage efficiencies. Guided by the H score (Supplementary Table 2), substituting the triplet repeat of the 7-amino-acid substrate with the corresponding 20-amino-acid natural sequence resulted in a 1.3- to 7.9-fold increase in cleavage-induced fluorescence (Fig. 4a). While GS linkers flanking the 7-amino-acid motif enhanced cleavage in some cases, the 20-amino-acid natural context consistently produced higher activity. This result suggests that H scores improve when we include the biochemical context surrounding the substrate sequence, which translates to improved proteolytic activity.

To determine whether ProSSpeC can be used to predict protease mutations that improve or preserve catalytic efficiency, we mutated five proteases, including the widely used TEVp and TVMVp, with similar H scores to their wildtype counterparts. We found that these mutants have comparable efficiencies to wildtype (Supplementary Fig. 4), demonstrating the ability of ProSSpeC to explore sequence space using H as a proxy for functionality.

We then tested whether ProSSpeC can predict how single amino-acid changes in the substrate sequence affect protease activity. We calculated changes in H scores (DeltaH) for point mutations of target substrates that predict improved (DeltaH < 0) or reduced (DeltaH > 0) proteolytic cleavage. We assayed nine single-site mutations in target substrates and found that DeltaH was a good predictor of cleavage outcomes (as measured by changes in fluorescence intensity, DeltaFluor) (Supplementary Fig. 5). For example, mutating the ScMVcs P3 position Phe→Tyr (DeltaH = -21.3) significantly improved the ScMVp-induced fluorescence in the assay (DeltaFluor = 0.54). By contrast, mutating the EAPVcs P2 Leu→His led to a large positive DeltaH (49.8) and significantly reduced EAPVp-induced fluorescence (DeltaFluor = -0.87). Furthermore, when comparing predicted increase or decrease of fluorescence across protease-substrate pairs, DeltaH is correlated with DeltaFluor (R = -0.69, p = 9.2 x 10-5) (Fig. 4b and Supplementary Table 3). These findings indicate that it is possible to model and predict protease activity against different substrates at single-residue resolution.

Encouraged by the resolving power of ProSSpeC, we tested whether it could prescribe functional proteases for a substrate that undergoes a single amino acid mutation. For example, ProSSpeC prescribed that EAPVp can cleave the SPV2cs when mutated at the P2 site (EPVYHQ^S→EPVYLQ^S; DeltaH = -49.8) (Fig. 4c, Supplementary Fig. 6, and Supplementary Table 3). We experimentally observed this crosstalk (DeltaFluor = 0.60), which was surprising given that mutSPV2cs (EPVYLQ^S) is markedly different from the native EAPVcs (ESVSLQ^T) that EAPVp cleaves (Fig. 2a). To examine the basis of this prediction, we analyzed the top protease–substrate residue couplings contributing to the lower H score. The ten strongest couplings all involved the mutated substrate P2 residue and numerous EAPVp residue positions, which were distributed across the protease (Fig. 4d and Supplementary Fig. 7). These mutagenesis experiments clearly demonstrate that ProSSpeC can resolve the effect of single-site substrate mutations on protease activity and identify crosstalk interactions between proteases.

ProSSpeC allowed us to identify functional proteases and engineer their crosstalk at single-residue resolution, which opens the door to unique biotechnological applications. To demonstrate this potential, we tested whether proteases can be used to selectively kill cells21,22 that harbor a point mutation in a target protein. We used and engineered Caspase3, an executioner caspase whose activation leads to cleavage of numerous endogenous proteins, ultimately leading to apoptosis. We tested whether EAPVp can synthetically trigger caspase-mediated apoptosis (i.e., synpoptosis) of mutant cells in a mixed population with three cell types: a wild-type cell (expressing Caspase3 with a wild-type SPV2cs: EPVYHQS); a mutant cell (expressing Caspase3 with mutSPV2cs: EPVYLQS); or a cell with no Caspase3 (as negative control). Flow cytometry of the mixed population showed that EAPVp can selectively trigger apoptosis in cells that expressed mutant Caspase3, but not in cells expressing wild-type or no Caspase3 (Fig. 5a; see Supplementary Fig. 8 for gating strategy). As expected, EAPVp efficiently cleaved and activated EAPVcs-containing Caspase3. Similarly, SPV2p only triggered apoptosis in cells with Caspase3 harboring SPV2cs (Fig. 5 b). This result demonstrates that engineered protease-substrate pairs can be used in biological circuits to detect and selectively trigger apoptosis of mutant cells in a mixed population.

Methods
Homolog search and alignment
Protease sequences were obtained using HMMER31. The hidden markov model profile (hmm seed) from the peptidase C4 family (InterPro PFAM accession: PF00863), which is the family that is used to describe NIa proteases, was used to run hmmsearch on the Uniprot databases (Sprot and TREMBL-February 2024 release). Sequences with greater than 45 contiguous gaps were removed from the alignment. The resulting protease alignment only contains sequences from Potyviridae. Substrate sequences were obtained by downloading the polyprotein sequences for the aligned and filtered proteases. To align the substrates such that the P1 position is at the same index, a manual alignment was created using the sequences denoted in Supplementary Table 4. The context for these sequences was extended to 20 amino acids based on previous findings26, which showed that surrounding amino acids could induce TEVp cleavage in non-canonical consensus sequences. The resulting hmm profile was then used to align the polyprotein sequences to the manual alignment. Substrate sequence composition is shown in Supplementary Fig. 10. Then, the proteases and substrates originating from the sam organism, as determined by Uniprot ID, were concatenated to each other to create the full alignment for specificity. This resulted in 3817 concatenated protease-substrate pairs, which are described in Supplementary Data 1. Another alignment was created for nonspecific signals by randomly scrambling the proteases and substrates. Both alignments were then used to calculate parameters for the ProSSpeC model.

Calculation and validation Protease-Substrate Specificity Calculator (ProSSpeC)
Direct coupling analysis (DCA) was performed using mfDCA9 on the canonical pairs of proteases and substrates or the scrambled pairs. To check the quality of the mfDCA results, we plotted the top contacts based on direct information (DI) value as shown in Supplementary Figs. 11 and 12. The learned DCA parameters from the canonical pairs were then used to compute the Potts Hamiltonian value:

Where H is the Hamiltonian value for canonical pairings, L is the length of the protease, L is the length of the substrate, e is the coupling between amino acids at positions i and j, and h is the local field for the amino acid at position i. Then, the scrambled pairs were used to calculate the Potts Hamiltonian value for the nonspecific signal:

The ProSSpeC model H score was then calculated by removing the nonspecific signal from the signal of the canonical pairs to yield the specific signal from the Hamiltonian:

In Fig. 2b, Fig. 3, Fig. 4a, Supplementary Fig. 1, and Supplementary Fig. 2, the substrate used in experimental assays was a triplet of the 7 amino acid consensus sequence. For these predictions, all substrate positions that are not P6-P1’ are masked and delineated using masked H :

H and fluorescence correlations and orthogonality metrics
Masked H values were clipped to the range [–180, –80], meaning that any score below –180 was set to –180 and any score above –80 was set to –80. This adjustment was made to make a fair comparison between fluorescence and H with the consideration that fluorescence intensity values cannot drop below 0. This cut-off is supported by the distribution of scrambled partners in comparison with natural partners. Masked H values were also multiplied by –1.0 to remove the inverse relationship, which is primarily important for the graph comparisons. In Fig. 3 correlation between masked H and fluorescence was calculated using Pearson R across all cognate and crosstalk pairs (R = 0.68; p = 1.40 × 10-31). A graph network similarity for protease-substrate specificity was also calculated using DeltaCon32. A bipartite graph was constructed such that each protease and each substrate was represented as a node with edges between each protease node and each substrate node, yielding 30 nodes and 225 edges. Edge weights were set to the normalized masked H values for the H graph and the normalized fluorescence values for the experimental graph. Similarity was calculated to be 0.676 with a Matusita distance of 0.479 after running with the following parameters: random_seed = 42, g = 20, ε = 0.01.

Phylogenetic Tree Construction
The WAG substitution model was used. Branch lengths, shown next to the branches (not drawn to scale), represent the number of substitutions per site. The initial tree for the heuristic search was selected by comparing neighbor-joining and maximum parsimony starting trees. Rate variation among sites was modeled using a discrete gamma distribution with four categories (+G, parameter = 1.8470), and 13.88% of sites were considered evolutionarily invariant (+I). Analyses were performed utilizing up to 12 parallel computing threads in MEGA1233.

Normalization of H orthogonality matrix
H scores were calculated for all pairwise protease–substrate combinations as previously described. To enable cross-protease comparison and reduce scale bias, raw H values were normalized relative to the corresponding diagonal (cognate-cleaving) interactions. For each element H , the normalized orthogonality score was computed by normalizing with respect to the cognate (diagonal) term H :

This yielded a complementary matrix—normalized (by division) (Supplementary Fig. 9)—representing relative specificity between non-cognate and cognate cleavage pairs. The matrix was visualized as heatmaps and subsequently used as input predictors for ROC analysis.

ROC analysis and performance metrics
Protease–substrate matrices (15 × 15) of predictors (H variants) and experimental readouts (fluorescence) were vectorized by row-major flattening to obtain paired lists of scores and labels. Ground-truth labels were defined from fluorescence using a prespecified threshold: positive (cleaving) if fluorescence ≥ 0.1. Because lower raw H values indicate stronger matches, raw H scores were sign-inverted prior to analysis (i.e., we used − H ); Normalized H were used as-is. Receiver operating characteristic (ROC) curves were computed in Python (v3.11) using scikit-learn (roc_curve, auc; v1.4). For each distinct score threshold, the true-positive rate (TPR) and false-positive rate (FPR) were calculated as:

where TP, true positive; FP, false positive; TN, true negative; FN, false negative. The area under the ROC curve (AUC) was obtained via trapezoidal integration (auc). An optimal decision threshold on the predictor was chosen by maximizing Youden’s J (TPR − FPR) across all thresholds.

Cleavage-induced fluorescence assay
We fused a cytoplasmic protein, arrestin beta-2, and sfCherry3C(1-10) (i.e., the first ten beta strands of the sfCherry3C fluorescent protein), such that cleavage of the latter allows it to bind with the eleventh beta strand sfCherry11M (fused to the histone H2B-sfGFP reporter) and reconstitute sCherry3C fluorescence 34. Measurements for protease-induced reconstitution of sCherry3C fluorescence were normalized to the reporter sfGFP fluorescence. To prevent previously known autoinhibition of Potyviral proteases17, we mutated serine in position 219 (in TEVp) to valine when present (see Supplementary Table 1 for homologous mutations in other proteases).

To benchmark the accuracy of the fluorescent protein reconstitution assay, the 31 pairs were also tested where each protease (and the nuclear reporter) was tagged to sfGFP and each substrate to mTagBFP2. Cells were then analyzed through flow cytometry (BD LSRFortessa™ Cell Analyzer). Cleavage-induced fluorescence intensities normalized against protease amount, substrate amount, or both yielded strong correlations (Pearson r = 0.8 (p = 1.8x10-26), 0.82 (p = 1.5x10-28), and 0.81 (p = 7.4x10-27), respectively) with the previous assay, suggesting that both assays reflect the actual efficiency of proteolysis (Supplementary Fig. 13).

Cloning
All plasmids use a pCS2+ backbone and were constructed using either Gibson assembly or site-directed mutagenesis (SDM) techniques (Supplementary Data 2). Protease sequences were synthesized (Integrated DNA Technologies, IDT), amplified using Phusion® High-Fidelity PCR Kit (New England Biolabs, NEB), and inserted into the pCS2+ backbone using In-Fusion® Snap Assembly Master Mix (Takara Bio). For point mutants, SDMs were performed using the Q5® Site-Directed Mutagenesis Kit (NEB). Oligonucleotide primers were also synthesized by IDT.

Cell culture
Human Embryonic Kidney (HEK) 293 cells (RRID:CVCL_0045) were a gift from the lab of Dr. Leonidas Bleris. Cells were cultured in a humidity-controlled incubator under standard culture conditions (5% CO , 37°C) using complete growth media consisting of Dulbecco's High Glucose Modified Eagles Medium (DMEM high glucose – Cytiva; catalog#SH30022FS), supplemented with 10% fetal bovine serum (FBS – Corning; catalog#MT35011CV), 1% penicillin—streptomycin (Corning; catalog#MT30002CI), and 1X MEM non-essential amino acid solution (NEAA – Sigma-Aldrich; catalog#M7145). Cells were passaged at 70-90% confluency. For subculturing, cells in T75 flasks were added with 3 mL of trypsin-ethylenediaminetetraacetic acid (EDTA) (Gibco; catalog#25200056) and incubated for 5 mins at 37°C. To inactivate trypsin, 5mL growth media was added. Cells were spun down (130xg, 5 mins), supernatant discarded, and cell pellet resuspended in 10mL of fresh growth media. A 1mL aliquot of the cell suspension was added for each T75 flask (Corning; catalog#07202000) containing 10mL of growth media for a 1:10 split ratio.

Transfection, imaging, and analysis
One day before transfection, HEK293 cells were seeded in 24-well plates at 75,000 cells/well. Transfections were done with the jetPRIME® reagent (Sartorius), following their recommended protocol. Each well was transfected with a ratio of 1:10:1 of protease:substrate:reporter (16 ng each of reporter and protease, 166 ng of substrate, and 50 ng of mock plasmid DNA not encoding any protein – adding up to 250 ng total DNA) mixed with 0.5 mL of jetPRIME® reagent in 25 mL of jetPRIME® buffer. All transfections were done this way except for Fig. 3b; where the same ratio of protease:substrate:reporter was followed except 300 ng of mock was added instead of 50 ng (adding up to 500 ng of total DNA), which were then mixed with 1 mL of jetPRIME® reagent in 50 mL of jetPRIME® buffer. Images were taken 24 hours after transfection using EVOS M5000 fluorescence microscope with a 20X objective (NA:0.45, WD:6.12 mm, Cat. no: AMEP4982). Fluorescence images were analyzed using FIJI. Fluorescence intensity values were calculated by subtracting background fluorescence from nuclei fluorescence and then normalizing output fluorescence (sfCherry3C) by dividing it with reporter fluorescence (sfGFP). For each condition, three biological replicates with two technical replicates each were done (corresponding to two wells). For each technical replicate, 30 nuclei were analyzed.

Flow Cytometry
Media from each well were discarded, and cells in each well were then detached by adding 250 μL of trypsin-EDTA and placed in the incubator (37°C, 5% CO ) for 5 minutes. Wells containing trypsin-detached cells were added with equal amount of fresh media to stop the reaction. Cells were transferred to 1.5mL tubes, spun down (130×g for 5 mins), washed once with ice-cold PBS, spun down again (130×g for 5 mins), and resuspended in 500 μL DMSEM high glucose + 2% FBS. Cells were transferred to 5-mL round bottom tubes through the cell strainer caps (352235, Falcon®) before running flow cytometry using the BD LSRFortessa™ Cell Analyzer. A minimum of 10,000 events were analyzed per condition.

Protease titration
Using the same transfection and analysis protocol, we tested five proteases across a range of plasmid concentrations (1.66 ng, 5 ng, 16.6 ng, and 50 ng) to evaluate their proteolytic efficiency. As protease expression increased, some proteases showed higher activity while others reached an early plateau (Supplementary Fig. 14). Notably, all proteases reached fluorescence saturation at 16.6 ng, the concentration used throughout this study. The resulting fluorescence therefore reflects apparent catalytic efficiencies that can be directly compared across proteases.

Structural modeling
We modeled the protease in complex with corresponding 7-amino-acid substrate using AlphaFold3 (alphafoldserver.com)15 to generate high-confidence structural predictions. The amino acid sequences of the protease and substrate were input as two different protein entities to capture protease-substrate interactions. The resulting top-ranked structural models were visualized using PyMOL (Version 3.1.3) to highlight the top ten couplings as well as the catalytic residue for each protease-substrate pair in Fig. 4d.

Selective synpoptosis of mutant cells
Caspase 3, whose natural cleavage site was replaced with either SPV2cs or mutSPV2cs, was co-transfected with either SPV2p or EAPVp into HEK293 cells. After 12 hours, cells were then stained with an apoptosis marker fluorophore-conjugated Annexin V, and the immunofluorescence intensities were measured using flow cytometry.

Synpoptosis transfection
One day before transfection, HEK293 cells were seeded in a 6-well plate at 300,000 cells/well. To generate a mixed population of cells transfected with a single Caspase3 type (EAPVp-cleavable Caspase3 tagged with mCherry; or SPV2p-cleavable Caspase3 tagged with sfGFP), liposome complexes were individually made for each plasmid. In Caspase3-only conditions, 20 ng of EAPVp-cleavable Caspase3 or SPV2p-cleavable Caspase3 was mixed with 230 ng of mock plasmid. Each mixture was then incubated with 1 μL of jetPRIME reagent in 50 μL of jetPRIME buffer. In conditions containing SPV2p and protease-cleavable Caspase3, 20 ng of EAPVp-cleavable Caspase3, 20 ng of SPV2p-cleavable Caspase3, or 60 ng of SPV2p-mTagBFP2 were separately mixed with mock plasmid to a total of 166ng for each two-component mixture. The same procedure was done with conditions containing EAPVp-mTagBFP2 (120 ng). Each mixture was then incubated with 0.66 μL of jetPRIME reagent in 33.3 μL of jetPRIME buffer. After a 10-min incubation, the complexes for each condition were added to each well. This procedure ensured that each cell in a mixed population expressed only either EAPVp-cleavable or SPV2p-cleavable Caspase3.

Synpoptosis staining and flow cytometry
Media from each well were individually collected into a 1.5-mL centrifuge tube. Cells in each well were detached using 1mL of trypsin-EDTA for 5 minutes. Trypsin-detached cells were transferred to previously collected media to stop the reaction. Cells were spun down (130×g for 5 mins), resuspended in fresh media, and placed in the incubator (37°C, 5% CO ) for 30 minutes to allow recovery. They were then spun down again, washed once with ice-cold PBS, and resuspended in Annexin V buffer (600 µL/well) with Annexin V stain (30 µL/well). Alexa Fluor 350-conjugated Annexin V (Cat. no. A23202, Invitrogen) was used for Caspase3 titration (see Supplementary Table 5) while Alexa Fluor 647-conjugated Annexin V (cat. no. A23204, Invitrogen) was used for the synpoptosis experiment (see Fig. 5 and Supplementary Fig. 8). Annexin V staining was done at 25°C for 15 minutes. Cells were then transferred to round bottom tubes through a cell strainer cap before running flow cytometry using the BD LSRFortessa™ Cell Analyzer. Gating shown in Supplementary Fig. 8. A total of 100,000 events were analyzed per condition.

Software
Protease and substrate sequences were gathered using Matlab 2023b and hmmer 3.4. In silico plots were generated using Python 3.11.5 (scipy 1.11.1, matplotlib 3.7.2, biopython 1.78, pandas 2.1.1, numpy 1.24.3, ete3 3.1.3). Fluorescence microscopy images were analyzed with ImageJ software (v 2.14.0/1.54f). For all flow cytometry experiments, flow cytometry data were processed using FlowJo software (v10.10.0). All figure schematics were generated or compiled using Adobe Illustrator (v. 29.6.1).

Statistics & Reproducibility
No statistical method was used to predetermine sample size. No data were excluded from the analyses. The experiments were not randomized. The Investigators were not blinded to allocation during experiments and outcome assessment.

All experiments were performed in biological triplicates, except for Fig. 3b as noted performed with n = 3-5 replicates. Where error bars are shown, results were also expressed as means ± SEM. Correlations were evaluated using Pearson’s correlation coefficient (r). Two-tailed (Fig. 3) or one-tailed (Fig. 4b, Supplementary Fig. 13) p-values were calculated, and correlations were considered statistically significant at p < 0.05.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "In this work, the authors present ProSSpeC (Protease Substrate Specificity Calculator), a coevolution-based statistical model derived from Direct Coupling Analysis (DCA) to predict the substrate specificity of Potyviridae NIa proteases. The authors built this model from over 3800 naturally paired protease-substrate sequences and used the Hamiltonian-based Hspec score to quantify specificity.\nThey validated the model through a fluorescent protein reconstitution assay in human cells, demonstrating a significant correlation between the model-predicted Hspec scores and observed proteolytic activity. Beyond validating native protease-substrate pairs, they applied the model to: Predict the impact of single-residue mutations on substrate cleavage; engineer orthogonal and crosstalk interactions; construct a synthetic protease-triggered cell death pathway (synpoptosis) that discriminates between mutant and wild-type cells.\nThis manuscript is rigorously constructed, starting from a clear guiding hypothesis, validated through extensive experimental assays, and expanded into multiple applications. It provides a conceptually solid framework that integrates bioinformatics and statistical mechanics, distinguishing itself from many black-box AI approaches in the field. The ProSSpeC model is elegant, interpretable, and practical, with promising implications for synthetic biology, therapeutic design, and enzyme engineering. I am supportive of the publication of this paper after addressing the following points:\na. Clarify the Hspec model in the Results section\nCurrently, the formulation and rationale of the Hspec model are only discussed in the Methods. A more detailed and accessible explanation of the model, including how the joint Hamiltonian is derived and why it works, should be included in the Results. This would allow readers to grasp the mechanistic intuition and significance of the model early on.\nb. Clarify the correlation analysis\nThe manuscript mentions a strong correlation between Hspec scores and cleavage-induced fluorescence (R = 0.79), but it's unclear how this correlation was calculated. Was every protease-substrate pair (including poorly active or non-cognate combinations) included in the analysis, or was the correlation assessed only on best-performing pairs? Please clarify the subset used and justify the approach.\nc. Explore mechanistic insight via structure-based modeling\nWhile the model performs impressively, the molecular mechanism underlying specificity remains abstract. I suggest that the authors select one representative protease-substrate pair (perhaps one with a well-characterized ΔHspec mutation) and perform structural modeling (e.g., docking or AlphaFold3-based complex prediction). This could shed light on how specific residue-residue interactions enhance or weaken specificity, thereby linking statistical couplings to physical mechanisms.\nd. Discuss out-of-cleavage-site effects\nIn the Discussion, the authors should comment on the limitation of focusing solely on the cleavage region. Mutations outside the P6–P1′ region (e.g., distal regions in the substrate protein) may alter local structure or dynamics, affecting cleavage. Acknowledging this limitation and suggesting future directions for extending the model (e.g., through structural context or experimental feedback loops) would strengthen the discussion.\ne. Address substrate promiscuity\nProSSpeC is designed around specificity prediction, but some proteases may exhibit promiscuous behavior—cleaving multiple substrates with moderate affinity. Could the authors comment on how the model handles such cases? Does a moderately negative Hspec score reliably indicate partial activity across multiple substrates?\nf. Clarify the CABMV case\nThe CABMV protease shows high cleavage efficiency toward its own substrate in the assay, but the Hspec score suggests potential promiscuity. This case deserves further discussion—does it suggest a decoupling between experimental signal strength and model specificity? Could this be a case of an energetically favorable but non-unique interaction? Understanding and explaining such discrepancies will help delineate the boundaries of the model.", "reviewer_2": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.", "reviewer_3": "In the manuscript “Identification and engineering of highly functional synthetic proteases in cells using co-evolutionary models,” the authors present Protease Substrate Specificity Calculator (ProSSpeC), a co-evolution-guided computational tool that predicts protease-substrate specificity. Using direct coupling analysis, ProSSpeC quantifies protease-substrate specificity via a Hamiltonian score. ProSSpeC was trained on over 3,000 natural NIa protease-substrate pairs from the Potyviridae family. It accurately predicted cleavage efficiency for many protease-substrate pairs, with Hamiltonian scores strongly correlated with experimental cleavage-induced fluorescence in human cells. Several proteases were found to outperform the commonly used TEV protease (TEVp) in human cells. Additionally, ProSSpeC predicted the effects of single-point mutations in the cleavage sites of substrates, enabling tuning of substrate specificity. The authors further used ProSSpeC to reprogram protease specificity and engineer crosstalk between a protease and a mutant non-cognate substrate. As a demonstration, they used a protease to selectively trigger apoptosis in cells expressing engineered caspase-3 containing the mutant cleavage site of another protease.\nOverall, this manuscript presents a co-evolutionary computational model supported by cell-based experimental validation. The model quantitatively predicts protease-substrate specificity at single-residue resolution. Integration of modeling with cell-based assays establishes a 'predict-validate' story. This study can offer programmable protease tools in synthetic biology, potentially enabling therapeutic applications such as selective cell killing.\nMajor Comments:\nIn the cellular fluorescence assay, it would be helpful if both the protease and the substrate were tagged with different fluorescent markers. That way, expression levels of both components can be monitored and controlled for by flow cytometry gating. Without rigorous quantification, it is difficult to tell whether differences in activity are really due to actual proteolysis efficiency or just variable expression or stability. Some data on protease and substrate levels would be necessary.\nIn the “ProSSpeC-guided engineering of nonnative cleavage specificity” section, it would add value to the paper if the authors systematically profile a library of substrate mutations for TEVp and TVMVp, which are widely used. Knowing how to improve cleavage efficiency for these proteases through targeted mutations would be valuable for many synthetic biology researchers. Right now, the focus is on ScMVp, SbMVp, and EAPVp, which are less common and could limit potential protease-based applications.\nProSSpeC offers quantitative Hamiltonian scores, but does not provide mechanistic insight or structural rationale for the most favorable protease-substrate pairs. Consider using AlphaFold-predicted complex structures or contact maps to rationalize top-scoring interactions.\nThe manuscript does not compare ProSSpeC against alternative protease-substrate specificity predictors, such as PROSPER, PeptideCutter, and iProt-Sub. A side-by-side comparison, even just for a subset of known pairs, would strengthen the claim that coevolutionary modeling adds value.\nThe observation that extending substrate sequences from 7 to 20 amino acids improves cleavage is interesting, but it is only shown for a few cases. Also, it would be good to clarify what the additional 13 amino acids are. Are they from natural sequence context, or are they synthetic linkers? Could flexible sequences like GGGS repeats also boost cleavage? Would the position of the 13 amino acids matter, such as on the N-terminus or the C-terminus of the 7-amino acid cleavage site?\nThe fluorescence signal for some protease-substrate pairs seems to plateau, potentially masking finer distinctions. It would be helpful to report absolute fluorescence ranges and dynamic range limits for the assay. Consider using titration experiments to better quantify proteolytic efficiency.\nThe ProSSpeC model is validated using a relatively limited number of 10×10 protease-substrate pairs. It would be great to see more pairs tested to better assess how generalizable the model is.\nMinor Comments:\nCan the authors comment on how this tool can or cannot be applied beyond the Potyviridae family? Would the same modeling approach generalize to other protease families?\nIt would be helpful to indicate where TEV protease lands on the distribution in Figure 2a—TEVp is a standard point of comparison, so readers will naturally look for it.\nMany experiments used only two independent replicates. N = 3 or more would be ideal.\nIn Figure 4b, it is not entirely clear what each dot represents. Clarifying that would help readers interpret the plot.\nFigure 4d could benefit from a few more control conditions where the cell expresses “EAPVp + Caspase3-EAPVcs” and “SPV2p + Caspase3-SPV2cs”.\nDefine synpoptosis more clearly on first mention.\nLastly, the authors wrote: “Our results highlight the potential of engineered proteases to interface with endogenous cellular pathways.” However, there is no demonstration of such interfacing in this study.", "reviewer_4": "Here the authors developed ProsspeC, a generative model that was trained using 3817 pairs of viral Nuclear Inclusion a (NIa) proteases and their cleavage peptide sequences. The model can generate a score (Hspec) for an input, with lower scores (negative scores) predicts higher chance of cleavage. The authors then validated this model using a cell-based assay based on reconstitution of fluorescence protein (successful cleavage releases a fragment, which enters the nucleus and complements the other fragment in the nucleus to form a complete fluorescent protein).\nThe writing is exceedingly short. This is not acceptable. The whole paper reads like a summary, instead of a full-length scientific paper. Lack of details renders it difficult to review this paper as many figures and authors’ interpretations are difficult to understand.\nThe study is limited to Nla protease family. No other types of proteases are examined. Authors need to tune down all broad claims of “proteases” and spell out clearly that the work is limited to Nla proteases in the title and abstract.\nThere is no “engineer protease” in this manuscript. Authors only tested mutated substrate sequences. The last paragraph of the “Discussion” needs to be corrected.\nIt is very strange that Fig. 1-4 are all covered by the first sub-title of the “Result” section, and the main figures only have these 4 figures. Many supplementary figures should be moved to main figures.\nThere are many cases that experimental data do not fit their Hspec prediction. Authors need to explain all these pairs that do not fit their models: Fig. 2b: SPV2p? PTVp? Fig. 3b: PPV versus KoMV substrate? Fig. S2: TEV: the score of -249 showed higher signal than -296. For PeMoV, -243 showed lower signal than -231? Fig. S3: for ScMV: the blue (-184) was predicted to be stronger than orange (-153), but the sfCherry3C signaling is higher for orange, how do you explain this? Fig. 4a: increase length to 20 did not drastically increase the degree of cleavage for ISMVp etc.\nFig. 3a: aren’t these pairs used to train the model? Why it is a surprise that they correlate?\nThe meaning of Fig. 4d is not clear: there is no engineering of protease. The authors just put a mutant substrate peptide sequence into caspase3. In this case, wouldn’t it be better to use the native substrate of this protease (EAPVp)?"}, "author_response": "Summary Response to Reviewers:\nWe thank the Reviewers for their constructive feedback on the manuscript. We performed substantial,\nadditional experiments, modeling, and clarifications on some conceptual and experimental details in the\ntext. We summarize the new results that further clarify or strengthen our previous claims:\n(1) In updated Fig. 3, we expanded the protease-substrate orthogonality matrix from 100 to 225\nexperimentally tested pairs. ProSSpeC H scores remained strongly correlated to the expanded\nspec\nexperimental dataset of on-target and crosstalk activity. Further, to rigorously quantify the predictive\ncapacity of H scores, we generated receiver operating characteristic (ROC) curves in new\nspec\nSupplementary Fig. S2. Raw or normalized H scores (AUC = 0.878; normalized AUC = 0.905) below a\nspec\nraw H cutoff of -84 can reliably discriminate between true cleaving vs non-cleaving pairs. We added\nspec\na new section in Methods to describe the ROC analysis.\n(2) In updated Fig. 4a, additional experiments on the substrate sequence context support the ProSSpeC\nprediction that the composition of the extended 20-a.a. substrate sequence (and not just linker length)\nsubstantially improves proteolytic activity for this family of proteases when compared to the minimal\n7-a.a. substrate motif.\n(3) In new Supplementary Fig. S3, we experimentally validated 12 additional protease triple-mutants\npredicted by the ProSSpeC model to have similar or more negative H scores than their wildtype\nspec\ncounterparts. These mutants exhibited cleavage efficiencies similar to wildtype. We note that mutations\ncan easily destroy protein function, but here we show that ProSSpeC-guided mutations efficiently\nmaintain the enzyme-substrate specificity relationship, allowing one to explore sequence space using\nthe H score as a predictor of functionality.\nspec\n(4) In new Fig. 4d, to offer possible structural insights on ProSSpeC predictions, we also performed\nAlphafold3 modeling to visualize how EAPVp effectively cleaves a non-canonical target (mutSPV2cs:\nSPV2cs with its P2 site mutated H à L). The top 50 amino-acid P2-couplings that contribute to a more\nfavorable H score for EAPVp-mutSPV2cs (compared to EAPVp-SPV2cs) are found in new\nspec\nSupplementary Fig. S6. We mapped the top 10 contributing residues of EAPVp onto Alphafold3 models\nof (non-cleaving pair) EAPVp-SPV2cs or (cleaving pair) EAPVp-mutSPV2cs complexes. A distributed\nnetwork of both active-site and epistatic residues contributes to the observed EAPVp-mutSPV2cs\ncrosstalk, which highlights the single-amino acid resolution and direct interpretability of ProSSpeC\nmodel predictions.\n(5) In new Supplementary Fig. S11, we benchmarked and validated our cleavage-induced fluorescence\nassay against new flow cytometry data from fluorescent protein-tagged protease and substrate, for all\n31 cognate pairs tested. Flow cytometry benchmarking further strengthened the robustness of the\nprevious fluorescence image analyses (Pearson r > 0.80). We thoroughly described this benchmarking\nassay in Methods.\n(6) In new Supplementary Fig. S12, we also performed protease titration experiments to demonstrate that\nall proteases reached fluorescence saturation at the protease amount used throughout this study, such\nthat readers can directly compare cleavage-induced fluorescence values between proteases.\n(7) All experiments were performed with independent replicates N = 3-5.\nWe have clarified the text where needed. We wrote specific responses to the Reviewers, as follows:\nSpecific Responses to Reviewers:\nReviewer #1 (Remarks to the Author):\nIn this work, the authors present ProSSpeC (Protease Substrate Specificity Calculator), a coevolution-based\nstatistical model derived from Direct Coupling Analysis (DCA) to predict the substrate specificity of Potyviridae\nNIa proteases. The authors built this model from over 3800 naturally paired protease-substrate sequences and\nused the Hamiltonian-based Hspec score to quantify specificity. They validated the model through a\nfluorescent protein reconstitution assay in human cells, demonstrating a significant correlation between the\nmodel-predicted Hspec scores and observed proteolytic activity. Beyond validating native protease-substrate\npairs, they applied the model to: Predict the impact of single-residue mutations on substrate cleavage; engineer\northogonal and crosstalk interactions; construct a synthetic protease-triggered cell death pathway\n(synpoptosis) that discriminates between mutant and wild-type cells. This manuscript is rigorously constructed,\nstarting from a clear guiding hypothesis, validated through extensive experimental assays, and expanded into\nmultiple applications. It provides a conceptually solid framework that integrates bioinformatics and statistical\nmechanics, distinguishing itself from many black-box AI approaches in the field. The ProSSpeC model is\nelegant, interpretable, and practical, with promising implications for synthetic biology, therapeutic design, and\nenzyme engineering.\nResponse 1.1:\nWe thank the Reviewer for supporting our work and appreciate the thoughtful comments and suggestions.\nWe hope your specific concerns are sufficiently addressed in this new revision, as discussed point-by-point\nbelow.\nReviewer #1:\nI am supportive of the publication of this paper after addressing the following points:\na. Clarify the Hspec model in the Results section\nCurrently, the formulation and rationale of the Hspec model are only discussed in the Methods. A more detailed\nand accessible explanation of the model, including how the joint Hamiltonian is derived and why it works,\nshould be included in the Results. This would allow readers to grasp the mechanistic intuition and significance\nof the model early on.\nResponse 1.2:\nDetails regarding the H model were added to the Results section, along with the rationale of its\nspec\nsignificance. We hope that the following changes provide readers with a stronger intuition of the model as\nthey read through the Results:\n“ProSSpeC is a coevolutionary model that leverages the covariation of NIa protease sequences along\nwith their aligned substrate sequences to calculate a Hamiltonian specificity score (H ), where the\nspec\nmore negative the H score is, the stronger the predicted specificity is. This specificity score is defined\nspec\nby the Potts Hamiltonian energy of canonical protease-substrate pairs (Equation 1), with the Potts\nHamiltonian energy of scrambled enzyme-substrate pairs (Equation 2) removed, to provide a score of\nsequence attributes unique to interacting pairs and not shared due to family-likeness or general\nsequence similarity (Equation 3). For both canonical and scrambled protease-substrate pairs, the Potts\nHamiltonian is defined as the sum of all parameters in a Boltzmann distribution, with couplings (e ) and\nij\nlocal fields (h) inferred using direct coupling analysis (DCA), specifically the mean-field formulation\ni\n(mfDCA)9. A similar rationale for using the DCA Hamiltonian to characterize molecular interactions has\nbeen applied to study specificity in two-component systems12, recognition in protein-RNA\ninteractions17, and to predict and engineer compatibility in hybrid transcription factors13,18. The\nProSSpeC model applies this approach to quantify interactions of the protease and the 7-amino acid\nsubstrate; we masked out substrate residues outside of this window in H calculations to yield the\nspec\nmasked H (Equation 4). The masked H was utilized for predictions where the designed experiment\nspec spec\ndid not include the extended substrate context. ProSSpeC allows us to quantitatively assess the\nspecificity of NIa proteases and their substrate motif through sequence alone, allowing us to filter and\nrank both natural NIa protease-substrate pairs as well as mutated ones prior to experimental testing.”\nReviewer #1:\nb. Clarify the correlation analysis\nThe manuscript mentions a strong correlation between Hspec scores and cleavage-induced fluorescence (R =\n0.79), but it's unclear how this correlation was calculated. Was every protease-substrate pair (including poorly\nactive or non-cognate combinations) included in the analysis, or was the correlation assessed only on best-\nperforming pairs? Please clarify the subset used and justify the approach.\nResponse 1.3:\nThe correlation calculation in this study took into account all the protease-substrate combinations that we\nhad tested; there was no pre-selection. In the Methods section, we have improved the language to assist\nthe reader in understanding the correlation calculation and added the following text to clarify our approach.\nAdditionally, in response to Reviewer 3, we now have more than doubled the number of protease-substrate\npairs that were experimentally tested (in new Fig. 3; 225 total pairs; Pearson R = 0.68).\nH and fluorescence correlations and orthogonality metrics\nspec\nMasked H values were clipped to the range [-180,-80], meaning that any score below –180 was\nspec\nset to –180 and any score above –80 was set to –80. This adjustment was made to make a fair\ncomparison between fluorescence and H with the consideration that fluorescence intensity\nspec\nvalues cannot drop below 0. This cut-off is supported by the distribution of scrambled partners in\ncomparison with natural partners. Masked H values were also multiplied by –1.0 to remove the\nspec\ninverse relationship, which is primarily important for the graph comparisons. In Fig. 3 correlation\nbetween masked H and fluorescence was calculated using PearsonR across all cognate and\nspec\ncrosstalk pairs (R = 0.68; p = 1.40 × 10-31). A graph network similarity for protease-substrate\nspecificity was also calculated using DeltaCon32 . A bipartite graph was constructed such that each\nprotease and each substrate was represented as a node with edges between each protease node\nand each substrate node, yielding 30 nodes and 225 edges. Edge weights were set to the\nnormalized masked H values for the H graph and the normalized fluorescence values for the\nspec spec\nexperimental graph. Similarity was calculated to be 0.676 with a Matusita distance of 0.479 after\nrunning with the following parameters: random_seed = 42, g = 20, ε = 0.01.\nReviewer #1:\nc. Explore mechanistic insight via structure-based modeling\nWhile the model performs impressively, the molecular mechanism underlying specificity remains abstract. I\nsuggest that the authors select one representative protease-substrate pair (perhaps one with a well-\ncharacterized ΔHspec mutation) and perform structural modeling (e.g., docking or AlphaFold3-based complex\nprediction). This could shed light on how specific residue-residue interactions enhance or weaken specificity,\nthereby linking statistical couplings to physical mechanisms.\nResponse 1.4:\nWe thank the Reviewer for their positive assessment of the performance of the ProSSpeC model. To further\nhighlight the interpretability of ProSSpeC at single-residue resolution, we have broken down the top 50\nresidue couplings (in new Supplementary Fig. S6) that contribute most to the model prediction that EAPVp-\nmutSPV2cs is a functionally cleaving pair, whereas EAPVp-SPV2cs is not (which was experimentally\ndemonstrated in Fig. 4c and new Fig. 5). As suggested by the Reviewer, we performed Alphafold3-based\nprediction of these two complexes to shed light on how statistical residue couplings might enhance or\nweaken specificity (new Fig. 4d).\n• Of the top 10 residue couplings for EAPVp-mutSPV2cs, three amino acids in the Alphafold3-\npredicted EAPVp structure are <10 Å from the mutated P2-site (H->L): W206, N173, and F179.\n• The top 1 coupling is a hydrophobic interaction between W206 and mutSPV2cs’ P2 L. The\nwildtype P2 is H, a positively charged residue which may be unfavorable for EAPVp’s W206.\n• The remaining top 7 couplings are from EAPVp residues > 10 Å away from P2, which imply that\ncleavage specificity also arises from epistatic interactions within a distributed network of\nresidues.\nWe have therefore added the following text in the Results section under the subsection “Engineering\ncrosstalk for selective synpoptosis”:\n“…To examine the basis of this prediction, we analyzed the top protease–substrate residue\ncouplings contributing to the lower H score. The ten strongest couplings all involved the\nspec\nmutated substrate P2 residue and numerous EAPVp residue positions, which were distributed\nacross the protease (Fig. 4d, and Supplementary Fig. S6).”\nAnd in the Discussion:\n“To demonstrate the interpretability of the model, we analyzed a case where SPV2cs, which is not\ncleaved by EAPVp, was predicted to become cleavable upon a single amino acid substitution (Fig.\n4c). We examined the top protease–substrate residue couplings contributing to this prediction. Of\nthe ten strongest couplings, only three EAPVp residues were within 10 Å of the mutated substrate\nposition (P2). The highest-ranked interaction was a hydrophobic contact between W206 and the\nmutant P2 Leu, whereas the wild-type P2 His, a positively charged residue, may be unfavorable for\nW206. The remaining seven top couplings involved residues located more than 10 Å from P2,\nsuggesting that cleavage specificity arises partly from long-range epistatic interactions within a\ndistributed residue network (Fig. 4d and Supplementary Fig. S6). These results illustrate how\nProSSpeC captures interpretable residue–residue dependencies underlying substrate recognition.”\nReviewer #1:\nd. Discuss out-of-cleavage-site effects\nIn the Discussion, the authors should comment on the limitation of focusing solely on the cleavage region.\nMutations outside the P6–P1′ region (e.g., distal regions in the substrate protein) may alter local structure or\ndynamics, affecting cleavage. Acknowledging this limitation and suggesting future directions for extending the\nmodel (e.g., through structural context or experimental feedback loops) would strengthen the discussion.\nResponse 1.5:\nWe thank the Reviewer for this suggestion, we now added more information, new analysis and updated Fig.\n4a, as well as commentary around the extended (20-a.a.) substrate context window in the Discussion:\n“We demonstrate that ProSSpeC, a quantitative model for predicting protease-substrate interaction,\nenables the identification and engineering of protease-substrate pairs with enhanced cleavage\nactivity in human cells. Unlike most previous analyses of NIa protease activity 23-25 , our ProSSpeC-\nguided experiments evaluated not only the 7-amino acid substrate but also its extended 20-amino\nacid context window, providing a more comprehensive assessment of protease-substrate specificity.\nPrevious work by Beaumont, et al.26 showed that extended substrate context matters for cleavage\nspecificity. Our family-wide analysis reveals that residues flanking the cleavage motif play a role in\ncleavage efficiency. Fig. 4a shows that the composition of amino acids (and not just length) outside\nof the P6-P1' region significantly improves catalytic efficiency (1.3 to 7.9-fold). This provides\nevidence that the extended context contributes to proteolytic cleavage. One potential explanation\nmay be that these regions have been subjected to evolutionary selection, potentially playing a role\nin the local structural dynamics of substrate recognition. Additionally, ProSSpeC can also predict\ncleavage activity at single-residue resolution (Fig. 4b,c), enabling programmable protease-\nsubstrate recognition. …”\nReviewer #1:\ne. Address substrate promiscuity\nProSSpeC is designed around specificity prediction, but some proteases may exhibit promiscuous behavior—\ncleaving multiple substrates with moderate affinity. Could the authors comment on how the model handles\nsuch cases? Does a moderately negative Hspec score reliably indicate partial activity across multiple\nsubstrates?\nf. Clarify the CABMV case\nThe CABMV protease shows high cleavage efficiency toward its own substrate in the assay, but the Hspec score\nsuggests potential promiscuity. This case deserves further discussion—does it suggest a decoupling between\nexperimental signal strength and model specificity? Could this be a case of an energetically favorable but non-\nunique interaction? Understanding and explaining such discrepancies will help delineate the boundaries of the\nmodel.\nResponse 1.6:\nIn the expanded orthogonality matrix of 225 protease-substrate combinations (updated Fig. 3a,b), ProSSpeC\nH scores correlated well with experimentally validated protease-substrate specificities (Pearson R = 0.68;\nspec\nDeltaCon similarity = 0.676). We can notice in the orthogonality matrix some off-diagonal regions that are\nsignals of crosstalk. The model does a good job finding these regions where a single protease can interact\nwith several peptides. To rigorously quantify this, we now generated a Receiver Operating Characteristic\n(ROC) curve for all 225 protease-substrate pairs in the orthogonality matrix. Using H score as the\nspec\npredictor and assuming an experiment threshold for successful cleavage at the normalized fluorescence of\n0.1, the area under the curve (AUC) = 0.87 and the optimal H cutoff = -84 (new Supplementary Fig.\nspec\nS2 and Methods). When H scores are far from the cutoff, ProSSpeC can reliably predict cleaving or\nspec\nnoncleaving pairs. However, H scores close to the cutoff are less reliable (most noticeably in CABMVp,\nspec\nEAPVp, and the newly added WMV2p and ACPp). Specifically, in these proteases, there are a lot of false\npositives being predicted. We noticed that there seem to be a protease sequence-dependent threshold, so\nwe normalized H score for each pair against the cognate H score (i.e., scores along the diagonal),\nspec spec\nby dividing with cognate score (compare similarity between new Supplementary Table S6 and updated Fig. 3b).\nThe AUC improved to 0.908, indicating that cognate-normalization improved discriminatory power (new\nSupplementary Fig. S2). The problem remains in that normalized H scores close to the cutoff are still less\nspec\nreliable. Nonetheless, normalization by cognate H values reduces the number of false positives for these\nspec\nfour proteases.\nOne possible explanation as to why there is a discrepancy between the predicted promiscuity and the\nexperimentally high substrate specificity for these proteases is that they may contain attributes that deviate\nfrom shared NIa protease family properties. ProSSpeC is trained to detect specific signals across all NIa\nprotease-substrate pairs, and while the model is reweighted to prevent bias from the training data, it is\nalways possible that these NIa protease-substrate pairs are not well represented (i.e., rare) within the\ndistributions found in natural sequence data. Expanding experimental coverage of low-homology\nsequences and retraining ProSSpeC accordingly could improve prediction accuracy and sequence diversity.\nThis is briefly discussed in the Discussion, but we have expanded to include this commentary about CABMVp.\nWe thus now have added this statement in the Results section under the new second subsection\n“Experimental validation of the ProSSpeC model”:\n“… Also, using these pairs and to quantify how the model captures crosstalk in off-diagonal\nelements , we generated a Receiver Operating Characteristic (ROC) curve with an Area Under the\nCurve (AUC) = 0.878 showing that the model exhibits strong discriminatory ability between true\ncleaving and non-cleaving pairs (H cutoff = -84).”\nspec\nAlso, these moderately negative H scores may denote that protease-substrate pairs do bind but do not\nspec\nresult in actual cleavage. We have added new text in Discussion describing this feature of the model and\nhow we can explain some observed crosstalk as well as correct some false positives. The new text reads:\n“One limitation shared by predictive models, including ProSSpeC, is the lower accuracy for\nsequences (or pairs of sequences) that significantly deviate from the natural distribution. For\nexample, ProSSpeC typically predicts low specificity for substrates with non-glutamine amino acids\nat the highly conserved P1 site. Another example is the predicted promiscuity but experimentally\nhigh substrate specificity for CABMVp in Fig. 3. These results suggest that CABMVp may contain\nattributes that deviate from the learned NIa protease family properties. Normalizing noncognate\nH scores relative to the cognate H score enhanced discriminatory performance and decreased\nspec spec\nthe incidence of false-positive predictions for CABMVp and other proteases. (Supplementary Fig.\nS2 and Supplementary Table S6). Still, to properly address this limitation, the low-homology\nsequence space needs to be experimentally explored. Retraining ProSSpeC on an expanded dataset\ncould enhance its predictive power and expand the sequence diversity of engineered protease-\nsubstrate pairs. Nonetheless, we have shown that the current model is effective at engineering pairs\nfor an expanded distribution based on natural sequence properties, supporting immediate use and\napplication of ProSSpeC.”\nReviewer #1 (Remarks on code availability):\nI failed to access the code as it shows \"this connection is not private\".\nResponse 1.7:\nThe website’s security certificate was temporarily expired. It has now been fixed, and a regularly scheduled\njob was added to the server to ensure that the security certificate is always valid.\nReviewer #2 (Remarks to the Author):\nI co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the\nNature Communications initiative to facilitate training in peer review and to provide appropriate recognition\nfor Early Career Researchers who co-review manuscripts.\nReviewer #2 (Remarks on code availability):\nClear and easy to use\nResponse 2.1:\nWe thank Reviewer #2 for their positive feedback.\nReviewer #3 (Remarks to the Author):\nIn the manuscript “Identification and engineering of highly functional synthetic proteases in cells using co-\nevolutionary models,” the authors present Protease Substrate Specificity Calculator (ProSSpeC), a co-evolution-\nguided computational tool that predicts protease-substrate specificity. Using direct coupling analysis,\nProSSpeC quantifies protease-substrate specificity via a Hamiltonian score. ProSSpeC was trained on over\n3,000 natural NIa protease-substrate pairs from the Potyviridae family. It accurately predicted cleavage\nefficiency for many protease-substrate pairs, with Hamiltonian scores strongly correlated with experimental\ncleavage-induced fluorescence in human cells. Several proteases were found to outperform the commonly\nused TEV protease (TEVp) in human cells. Additionally, ProSSpeC predicted the effects of single-point\nmutations in the cleavage sites of substrates, enabling tuning of substrate specificity. The authors further used\nProSSpeC to reprogram protease specificity and engineer crosstalk between a protease and a mutant non-\ncognate substrate. As a demonstration, they used a protease to selectively trigger apoptosis in cells expressing\nengineered caspase-3 containing the mutant cleavage site of another protease.\nOverall, this manuscript presents a co-evolutionary computational model supported by cell-based\nexperimental validation. The model quantitatively predicts protease-substrate specificity at single-residue\nresolution. Integration of modeling with cell-based assays establishes a 'predict-validate' story. This study can\noffer programmable protease tools in synthetic biology, potentially enabling therapeutic applications such as\nselective cell killing.\nResponse 3.1:\nWe thank the reviewer for taking the time to review our work. We are indeed excited to share the\nprogrammable proteases to the broader synthetic biology community. We hope we have sufficiently\naddressed your specific concerns in this new revision, as discussed point-by-point below.\nReviewer #3\nMajor Comments:\nIn the cellular fluorescence assay, it would be helpful if both the protease and the substrate were tagged with\ndifferent fluorescent markers. That way, expression levels of both components can be monitored and controlled\nfor by flow cytometry gating. Without rigorous quantification, it is difficult to tell whether differences in activity\nare really due to actual proteolysis efficiency or just variable expression or stability. Some data on protease\nand substrate levels would be necessary.\nResponse 3.2:\nWe have now done substantial, additional experiments wherein each protease (and reporter) was tagged\nwith sfGFP and each substrate was tagged with mTagBFP2. New Supplementary Fig. S11 shows the\ncorrelation between cleavage-induced fluorescence quantified through fluorescence microscopy imaging\n(Method 1: without fluorescent protein tags for either protease or substrate) and cleavage-induced\nfluorescence quantified through flow cytometry (Method 2: fluorescent protein-tagged protease and\nsubstrate). We observed strong correlations between both Methods, whether cleavage-induced\nfluorescence was normalized to protease amount (Pearson r = 0.80), to substrate amount (Pearson r = 0.82),\nor to both (Pearson r = 0.81). We have now added these in the Methods section under “Cleavage-induced\nfluorescence assay”:\n“To benchmark the accuracy of the fluorescent protein reconstitution assay, the 31 pairs were also\ntested where each protease (and the nuclear reporter) was tagged to sfGFP and each substrate to\nmTagBFP2. Cells were then analyzed through flow cytometry (BD LSRFortessa™ Cell Analyzer).\nCleavage-induced fluorescence intensities normalized against protease amount, substrate amount,\nor both yielded strong correlations (Pearson r = 0.8 (p = 1.9x10-26), 0.82 (p = 1.5x10-28), and 0.81 (p\n= 7.8x10-27), respectively) with the previous assay, suggesting that both assays reflect the actual\ne actual\nefficiency of proteolysis (Supplementary Fig. S11).”\nReviewer #3\nIn the “ProSSpeC-guided engineering of nonnative cleavage specificity” section, it would add value to the paper\nif the authors systematically profile a library of substrate mutations for TEVp and TVMVp, which are widely\nused. Knowing how to improve cleavage efficiency for these proteases through targeted mutations would be\nvaluable for many synthetic biology researchers. Right now, the focus is on ScMVp, SbMVp, and EAPVp, which\nare less common and could limit potential protease-based applications.\nResponse 3.3:\nWe thank the reviewer for this suggestion. We now have tested 12 additional triple mutants for TEVp,\nTVMVp, as well as the lesser known LMVp, TVBMVp, PMAVp, and WMV2p (with similar or better H\nspec\nagainst their wild-type substrates). Supplementary Fig. S3 show that these mutants have essentially similar\ncleavage efficiencies compared to their wildtype counterparts. We note that mutations can easily destroy\nprotein function, but here we show that ProSSpeC-guided mutations efficiently maintain the enzyme-\nsubstrate specificity relationship, allowing one to explore sequence space using the H score as a\nspec\nprediction of functionality. The Hamiltonian score is able to quantify from a larger potential sequence space:\nthose mutations that maintain amino acid couplings even if the sequence is different. Evolutionarily, this\nrepresents a space of neutral or nearly-neutral mutations that is required to provide robustness and\nmaintain functionality even when mutations occur. We have added the following statement under the\n“ProSSpeC-guided engineering of nonnative cleavage specificity” section of the main text:\n“To determine whether ProSSpeC can be used to predict protease mutations that improve or\npreserve catalytic efficiency, we mutated five proteases, including the widely used TEVp and TVMVp,\nwith similar H scores to their wildtype counterparts. We found that these mutants have\nspec\ncomparable efficiencies to wildtype (Supplementary Fig. S3), demonstrating the ability of\nProSSpeC to explore sequence space using H as a proxy for functionality.“\nspec\nReviewer #3:\nProSSpeC offers quantitative Hamiltonian scores, but does not provide mechanistic insight or structural\nrationale for the most favorable protease-substrate pairs. Consider using AlphaFold-predicted complex\nstructures or contact maps to rationalize top-scoring interactions.\nResponse 3.4:\nSimilar to Response 1.4 and as also suggested by this Reviewer, we performed Alphafold3-based prediction\nof these two complexes to shed light on how statistical residue couplings might enhance or weaken\nspecificity (new Fig. 4d).\n• Of the top 10 residue couplings for EAPVp-mutSPV2cs, three amino acids in the Alphafold3-\npredicted EAPVp structure are <10 Å from the mutated P2-site (H->L): W206, N173, and F179.\n• The top 1 coupling is a hydrophobic interaction between W206 and mutSPV2cs’ P2 L. The\nwildtype P2 is H, a positively charged residue which may be unfavorable for EAPVp’s W206.\n• The remaining top 7 couplings are from EAPVp residues > 10 Å away from P2, which imply that\ncleavage specificity also arises from epistatic interactions within a distributed network of\nresidues.\nWe have therefore added the following text in the Results section under the subsection “Engineering\ncrosstalk for selective synpoptosis”:\n“…To examine the basis of this prediction, we analyzed the top protease–substrate residue\ncouplings contributing to the lower H score. The ten strongest couplings all involved the\nspec\nmutated substrate P2 residue but different EAPVp positions, which were distributed across the\nprotease (Fig. 4d, and Supplementary Fig. S6).”\nAnd in the Discussion:\n“To demonstrate the interpretability of the model, we analyzed a case where SPV2cs, which is not\ncleaved by EAPVp, was predicted to become cleavable upon a single amino acid substitution (Fig.\n4c). We examined the top protease–substrate residue couplings contributing to this prediction. Of\nthe ten strongest couplings, only three EAPVp residues were within 10 Å of the mutated substrate\nposition (P2). The highest-ranked interaction was a hydrophobic contact between W206 and the\nmutant P2 Leu, whereas the wild-type P2 His, a positively charged residue, may be unfavorable for\nW206. The remaining seven top couplings involved residues located more than 10 Å from P2,\nsuggesting that cleavage specificity arises partly from long-range epistatic interactions within a\ndistributed residue network (Fig. 4d and Supplementary Fig. S6). These results illustrate how\nProSSpeC captures interpretable residue–residue dependencies underlying substrate recognition.”\nReviewer #3:\nThe manuscript does not compare ProSSpeC against alternative protease-substrate specificity predictors, such\nas PROSPER, PeptideCutter, and iProt-Sub. A side-by-side comparison, even just for a subset of known pairs,\nwould strengthen the claim that coevolutionary modeling adds value.\nResponse 3.5:\nWe thank the Reviewer for this thoughtful comment and recognize that several computational models exist\nfor general protease cleavage, however we found them to not be applicable to our specific evaluation of\nNIa proteases:\n1. PeptideCutter only evaluates the Tobacco Etch Virus protease (TEVp) from the NIa peptidase family and\ndoes not give any indication of specificity, other than TEVp will cut at the glutamine residue if present\nat P1 and not otherwise. Mutating one residue in the 7-amino acid substrate motif sequence has no\neffect on this prediction, mutating more than one amino acid automatically changes the prediction to\nno cutting. Furthermore, mutations in the extended context have no impact on these TEVp predictions.\n2. ProsperousPlus (now the location of PROSPER and iProtSub) also hosts many protease predictions, but\nit does not support the NIa protease family, including key enzymes like TEVp.\nNonetheless, we do recognize the work of other scientists in this field, so we have chosen to highlight other\nuseful resources across protease families in our Discussion:\n“While other models for protease cleavage exist, they focus on protease families other than NIa\nproteases27-30, and more recently, only on TEVp25. By contrast, our ProSSpeC model provides full\nsupport for the entire NIa protease family (including TEVp) and leverages coevolutionary\ninformation to guide NIa protease-substrate specificity predictions, allowing for a comprehensive\nanalysis of the underlying sequence space of both enzyme and substrate. We believe that our tool\nfills an important research niche that complements existing resources, providing computational\nsupport that did not previously exist.”\nReviewer #3:\nThe observation that extending substrate sequences from 7 to 20 amino acids improves cleavage is interesting,\nbut it is only shown for a few cases. Also, it would be good to clarify what the additional 13 amino acids are.\nAre they from natural sequence context, or are they synthetic linkers?\nResponse 3.6:\nThe amino acids for the extended context (20 amino acids) comes from the natural sequence of the Potyviral\npolyprotein that surrounds the 7-amino acid motif. Supplementary Table S7 contains all the extended\nsubstrate sequences for all experiments. We have updated the legend in updated Fig. 4a to clarify the exact\nsubstrate sequences used in the different test conditions.\nReviewer #3:\nCould flexible sequences like GGGS repeats also boost cleavage? Would the position of the 13 amino acids\nmatter, such as on the N-terminus or the C-terminus of the 7-amino acid cleavage site?\nResponse 3.7:\nWe thank the reviewer for this comment. We have now added new data on protease activity against GS\nlinker-flanked cleavage motifs in updated Fig. 4a. While there are cases where this shows improvement over\nthe 3x7 amino acid cleavage motifs, most cases show that the extended motif (20 amino acid natural\ncleavage sequence) does improve cleavage efficiency. We have added the following statement in the Results\nsection under the new second subsection “Experimental validation of the ProSSpeC model”:\n“….Guided by the H score (Supplementary Table S4), substituting the triplet repeat of the 7-\nspec\namino-acid substrate with the corresponding 20-amino-acid natural sequence resulted in a 1.3- to\n7.9-fold increase in cleavage-induced fluorescence (Fig. 4a). While GS linkers flanking the 7-amino-\nacid motif enhanced cleavage in some cases, the 20-amino-acid natural context consistently\nproduced higher activity….”\nReviewer #3:\nThe fluorescence signal for some protease-substrate pairs seems to plateau, potentially masking finer\ndistinctions. It would be helpful to report absolute fluorescence ranges and dynamic range limits for the assay.\nConsider using titration experiments to better quantify proteolytic efficiency.\nResponse 3.8:\nSupplementary Fig. S1 shows quite a wide dynamic range such that the relative efficiencies are sufficiently\ncaptured by the assay. In new Supplementary Fig. S11, we now present additional cleavage data on all 31\nproteases and show that the cleavage-induced fluorescence values are well below the maximum possible\nsfCherry/sfGFP achieved by the positive control condition. In new Supplementary Fig. S12, we have also\ndone dose titration on 5 representative proteases, which show that the protease amount we used in all\nother experiments (1X = 16.6 ng) is indeed at saturation, indicating that the cleavage-induced fluorescence\nvalues reflect their relative activities. We have added this to the Methods section under the subsection\n“Protease titration”:\n“Using the same transfection and analysis protocol, we tested five proteases across a range of\nplasmid concentrations (1.66 ng, 5 ng, 16.6 ng, and 50 ng) to evaluate their proteolytic efficiency.\nAs protease expression increased, some proteases showed higher activity while others reached an\nearly plateau (Supplementary Fig. S12). Notably, all proteases reached fluorescence saturation at\n16.6 ng, the concentration used throughout this study. The resulting fluorescence therefore reflects\napparent catalytic efficiencies that can be directly compared across proteases.”\nReviewer #3:\nThe ProSSpeC model is validated using a relatively limited number of 10×10 protease-substrate pairs. It would\nbe great to see more pairs tested to better assess how generalizable the model is.\nResponse 3.9:\nWe thank the reviewer for this comment. We have more than double the pairs tested from 100 to 225\nprotease-substrate pairs, giving us a 15x15 orthogonality matrix (in updated Fig. 3).\nReviewer #3:\nMinor Comments:\nCan the authors comment on how this tool can or cannot be applied beyond the Potyviridae family?\nWould the same modeling approach generalize to other protease families?\nResponse 3.10:\nThis tool only applies to the Potyviridae family. The method can certainly be applied to other protease\nfamilies if retrained on sequence datasets for other families. The ProSSpeC model itself cannot generalize\nacross protease families without retraining.\nReviewer #3:\nIt would be helpful to indicate where TEV protease lands on the distribution in Fig. 2a—TEVp is a standard\npoint of comparison, so readers will naturally look for it.\nResponse 3.11:\nWe have updated Fig. 2a to indicate where TEVp is.\nReviewer #3:\nMany experiments used only two independent replicates. N = 3 or more would be ideal.\nResponse 3.12:\nAll experiments are now N >= 3.\nReviewer #3:\nIn Fig. 4b, it is not entirely clear what each dot represents. Clarifying that would help readers interpret the plot.\nResponse 3.13:\nWe have clarified what each dot represents in the caption for Fig. 4b:\n“b, DeltaH (ΔH ) versus DeltaFluor (Δ fluorescence) for each dot, which represents the difference\nspec\nbetween a protease cleaving one substrate versus the same protease cleaving another substrate\n(see also Supplementary Table S5).”\nReviewer #3:\nFig. 4d could benefit from a few more control conditions where the cell expresses “EAPVp + Caspase3-EAPVcs”\nand “SPV2p + Caspase3-SPV2cs”.\nResponse 3.14:\nWe have placed all synpoptosis experiment results in new Fig. 5 including positive controls “EAPVp +\nCaspase3-EAPVcs” and “SPV2p + Caspase3-SPV2cs”.\nReviewer #3:\nDefine synpoptosis more clearly on first mention.\nResponse 3.15:\nWe now have clarified the definition of synpoptosis under the subsection “Engineering crosstalk for selective\nsynpoptosis”:\n“…We tested whether EAPVp can synthetically trigger caspase-mediated apoptosis (i.e.,\nsynpoptosis)…”\nReviewer #3:\nLastly, the authors wrote: “Our results highlight the potential of engineered proteases to interface with\nendogenous cellular pathways.” However, there is no demonstration of such interfacing in this study.\nResponse 3.16:\nWe apologize for the confusion. We now added a statement clarifying that the engineered Caspase3 used\nin the study interfaces with the endogenous cellular pathway that result in apoptosis (2nd paragraph of\n“Engineering crosstalk for selective synpoptosis” subsection:\n“…To demonstrate this potential, we tested whether proteases can be used to selectively kill cells21,22\nthat harbor a point mutation in a target protein. We used and engineered Caspase3, an executioner\ncaspase whose activation leads to cleavage of numerous endogenous proteins, ultimately leading\nto apoptosis.”\nReviewer #3 (Remarks on code availability):\nn/a\nReviewer #4 (Remarks to the Author):\nHere the authors developed ProsspeC, a generative model that was trained using 3817 pairs of viral Nuclear\nInclusion a (NIa) proteases and their cleavage peptide sequences. The model can generate a score (Hspec) for\nan input, with lower scores (negative scores) predicts higher chance of cleavage. The authors then validated\nthis model using a cell-based assay based on reconstitution of fluorescence protein (successful cleavage\nreleases a fragment, which enters the nucleus and complements the other fragment in the nucleus to form a\ncomplete fluorescent protein).\nThe writing is exceedingly short. This is not acceptable. The whole paper reads like a summary, instead of a\nfull-length scientific paper. Lack of details renders it difficult to review this paper as many Figures and authors’\ninterpretations are difficult to understand.\nResponse 4.1:\nWe have now included additional details to better explain concepts and Figures in the paper. Details\nregarding the H model were added to the Results section, along with the rationale of its\nspec\nsignificance. We hope that the following changes provide readers with a stronger intuition of the model as\nthey read through the Results:\n“ProSSpeC is a coevolutionary model that leverages the covariation of NIa protease sequences along\nwith their aligned substrate sequences to calculate a Hamiltonian specificity score (H ), where the\nspec\nmore negative the H score is, the stronger the predicted specificity is. This specificity score is defined\nspec\nby the Potts Hamiltonian energy of canonical protease-substrate pairs (Equation 1), with the Potts\nHamiltonian energy of scrambled enzyme-substrate pairs (Equation 2) removed, to provide a score of\nsequence attributes unique to interacting pairs and not shared due to family-likeness or general\nsequence similarity (Equation 3). For both canonical and scrambled protease-substrate pairs, the Potts\nHamiltonian is defined as the sum of all parameters in a Boltzmann distribution, with couplings (e ) and\nij\nlocal fields (h) inferred using direct coupling analysis (DCA), specifically the mean-field formulation\ni\n(mfDCA)9. A similar rationale for using the DCA Hamiltonian to characterize molecular interactions has\nbeen applied to study specificity in two-component systems12, recognition in protein-RNA\ninteractions17, and to predict and engineer compatibility in hybrid transcription factors13,18. The\nProSSpeC model applies this approach to quantify interactions of the protease and the 7-amino acid\nsubstrate; we masked out substrate residues outside of this window in H calculations to yield the\nspec\nmasked H (Equation 4). The masked H was utilized for predictions where the designed experiment\nspec spec\ndid not include the extended substrate context. ProSSpeC allows us to quantitatively assess the\nspecificity of NIa proteases and their substrate motif through sequence alone, allowing us to filter and\nrank both natural NIa protease-substrate pairs as well as mutated ones prior to experimental testing.”\nReviewer #4:\nThe study is limited to Nla protease family. No other types of proteases are examined. Authors need to tune\ndown all broad claims of “proteases” and spell out clearly that the work is limited to Nla proteases in the title\nand abstract.\nResponse 4.2:\nWe have gone through the text and specified that we are looking at NIa protease-substrate interactions in\nseveral places where clarification was needed.\nReviewer #4:\nThere is no “engineer protease” in this manuscript. Authors only tested mutated substrate sequences. The last\nparagraph of the “Discussion” needs to be corrected.\nResponse 4.3:\nWe thank the reviewer for this comment and recognize that we needed to highlight more information on\nthe proteases themselves. When we performed in silico mutation of the proteases, we found very few\nmutants that had H scores that are similar to, if not better than, their wild-type counterparts. Of those\nspec\nthat we did find and test, they still maintained similar functionality to their wild-type counterparts. This is\nnow shown in new Supplementary Fig. S3, where five proteases were mutated at three sites, and we saw\nlittle to no difference in cleavage-induced fluorescence. This work will be expanded in the future, where we\nwill attempt to mutate both enzyme and substrate in conjunction with each other to access new sequence\nspaces. Thus, we have now added this paragraph in the Results section under the subsection “ProSSpeC-\nguided engineering of nonnative cleavage specificity”:\n“To determine whether ProSSpeC can be used to predict protease mutations that improve or\npreserve catalytic efficiency, we mutated five proteases, including the widely used TEVp and TVMVp,\nwith similar H scores to their wildtype counterparts. We found that these mutants have\nspec\ncomparable efficiencies to wildtype (Supplementary Fig. S3), demonstrating the ability of\nProSSpeC to explore sequence space using H as a proxy for functionality.”\nspec\nReviewer #4:\nIt is very strange that Fig. 1-4 are all covered by the first sub-title of the “Result” section, and the main Figures\nonly have these 4 Figures. Many supplementary Figures should be moved to main Figures.\nResponse 4.4:\nWe have broken up Fig. 4 to create new Figs. 4 and 5. In new Fig. 4d, we included Alphafold3 models that\nattempt to explain why the engineered EAPVp-mutSPVcs pair works. In new Fig. 5, we also added additional\nexperimental protease-substrate combinations to comprehensively showcase the synpoptosis use case, as\nwell as to improve the flow of the text.\nReviewer #4:\nThere are many cases that experimental data do not fit their Hspec prediction. Authors need to explain all\nthese pairs that do not fit their models: Fig. 2b: SPV2p? PTVp? Fig. 3b: PPV versus KoMV substrate? Fig. S2:\nTEV: the score of -249 showed higher signal than -296. For PeMoV, -243 showed lower signal than -231? Fig.\nS3: for ScMV: the blue (-184) was predicted to be stronger than orange (-153), but the sfCherry3C signaling is\nhigher for orange, how do you explain this? Fig. 4a: increase length to 20 did not drastically increase the degree\nof cleavage for ISMVp etc.\nResponse 4.5:\nWhile the model is not perfect, H scores and experimental fluorescence values are shown to be\nspec\ncorrelated (Fig. 3, Pearson R = 0.68; Fig. 4b, Pearson R = -0.67). Discrepancies between predicted and\nexperimental values may be attributed to the model not being able to delineate between sequence effects\non substrate binding vs actual substrate cleavage (e.g., a more negative H score may indicate higher\nspec\nbinding affinity but not necessarily increased cleavage efficiency). Also, small differences in ProSSpeC scores\nare less reliable than large differences. Additionally, due to the sequence-based nature of the model, NIa\nenzyme-substrate pairs which do not adhere to the natural distributions found across the family are less\nreliable. This is addressed in the Discussion and is already a point of interest for future work, where we can\nimprove our modeling capabilities by expanding sequence pair novelty in the training data.\n“One limitation shared by predictive models, including ProSSpeC, is the lower accuracy for\nsequences (or pairs of sequences) that significantly deviate from the natural distribution. For\nexample, ProSSpeC typically predicts low specificity for substrates with non-glutamine amino acids\nat the highly conserved P1 site. Another example is the predicted promiscuity but experimentally\nhigh substrate specificity for CABMVp in Fig. 3. These results suggest that CABMVp may contain\nattributes that deviate from the learned NIa protease family properties. Normalizing noncognate\nH scores relative to the cognate H score enhanced discriminatory performance and decreased\nspec spec\nthe incidence of false-positive predictions for CABMVp and other proteases. (Supplementary Fig.\nS2 and Supplementary Table S6). Still, to properly address this limitation, the low-homology\nsequence space needs to be experimentally explored. Retraining ProSSpeC on an expanded dataset\ncould enhance its predictive power and expand the sequence diversity of engineered protease-\nsubstrate pairs. Nonetheless, we have shown that the current model is effective at engineering pairs\nfor an expanded distribution based on natural sequence properties, supporting immediate use and\napplication of ProSSpeC.”\nWe now also generated a Receiver Operating Characteristic (ROC) curve for all 225 protease-substrate pairs\ntested. Using H score as the predictor and assuming an experiment threshold for successful cleavage at\nspec\nthe normalized fluorescence of 0.1, the area under the curve (AUC) = 0.87 and the optimal H cutoff = -\nspec\n84 (new Supplementary Fig. S2 and Methods). We thus now have added this statement in the Results\nsection under the new second subsection “Experimental validation of the ProSSpeC model”:\n“… Also, using these pairs and to quantify how the model captures crosstalk in off-diagonal\nelements, we generated a Receiver Operating Characteristic (ROC) curve with an Area Under the\nCurve (AUC) = 0.878 showing that the model exhibits strong discriminatory ability between true\ncleaving and non-cleaving pairs (H cutoff = -84).”\nspec\nReviewer #4:\nFig. 3a: aren’t these pairs used to train the model? Why it is a surprise that they correlate?\nResponse 4.6:\nProSSpeC is an unsupervised model trained on natural sequences. The correlation is between the\nProSSpeC's predicted score and the experimental functional cleavage assay. Beyond the cognate protease-\nsubstrate pairs, the unsupervised model is also able to capture (off-diagonal) crosstalk interactions.\nReviewer #4:\nThe meaning of Fig. 4d is not clear: there is no engineering of protease. The authors just put a mutant substrate\npeptide sequence into caspase3. In this case, wouldn’t it be better to use the native substrate of this protease\n(EAPVp)?\nResponse 4.7:\nWe apologize for any confusion in Fig. 4d (now new Fig. 5) and have added to the text that Caspase3 is the\nprotease, which we have engineered in this specific plot. This figure conceptualizes whether or not point\nmutations in a target protein can be detected by a designer protease in a mixed population of cells. The\nexperiments represent a proof-of-concept for our ultimate goal of creating mutation-detecting proteases.\nWe hope to leverage the findings from this work to expand both protease sequences and substrate\nsequences into a library of targetable sequence motifs and engineered proteases, while also understanding\noff-target effects. Here, we were able to leverage off-target effects for synpoptosis.\nREVIEWERS' COMMENTS\nReviewer #1 (Remarks to the Author):\nThe authors have addressed all my critiques in the revised manuscript. I am supportive of publication.\nReviewer #3 (Remarks to the Author):\nThe revision addresses all of our concerns.\nRESPONSE TO REVIEWERS:\nWe thank the Reviewers for their positive feedback on the revised manuscript."}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors have addressed all my critiques in the revised manuscript. I am supportive of publication.", "reviewer_3": "The revision addresses all of our concerns."}, "author_response": "We thank the Reviewers for their positive feedback on the revised manuscript."}]

Decision letter