Physics-based generative model of curvature sensing peptides; distinguishing sensors from binders

Proteins can specifically bind to curved membranes through curvature-induced hydrophobic lipid packing defects. The chemical diversity among such curvature “sensors” challenges our understanding of how they differ from general membrane “binders” that bind without curvature selectivity. Here, we combine an evolutionary algorithm with coarse-grained molecular dynamics simulations (Evo-MD) to resolve the peptide sequences that optimally recognize the curvature of lipid membranes. We subsequently demonstrate how a synergy between Evo-MD and a neural network (NN) can enhance the identification and discovery of curvature sensing peptides and proteins. To this aim, we benchmark a physics-trained NN model against experimental data and show that we can correctly identify known sensors and binders. We illustrate that sensing and binding are phenomena that lie on the same thermodynamic continuum, with only subtle but explainable differences in membrane binding free energy, consistent with the serendipitous discovery of sensors.


INTRODUCTION
The recognition of curved regions of lipid bilayer membranes by proteins plays a key role in many biological processes, such as vesicular transport, fusion, and fission (1,2). This preferred binding to curved membranes is called curvature sensing and is driven by the outer leaflet of the curved bilayer membrane being stretched, which causes defects in the packing of the polar lipid head groups. Apolar amino acids of proteins can complement the now exposed hydrophobic tails within these lipid packing defects, negating their energetic penalty and resulting in a thermodynamic driving force (Fig. 1A).
Besides fundamental biological importance, curvature selectivity has been proposed as a potential avenue for the development of broad-spectrum antiviral peptides that leverage the difference in curvature between the membranes of small enveloped viruses and the essentially flat host cell membrane (3)(4)(5)(6). However, the extremely serendipitous discovery and resulting rarity of curvature selective peptides obstruct the utilization of state-of-the-art data sciencedriven generative models, like recent work on the discovery of antimicrobial peptides (7). Consequently, an efficient computational strategy for accelerating the discovery of curvature sensing peptides is still lacking.
Many natural curvature sensing proteins feature an amphipathic helix (AH). AHs have a polar face that interacts with the solvent and the lipid head groups and an apolar face that interacts with the hydrophobic lipid tails. Beyond this shared structural amphipathicity, the chemical composition of AHs is highly diverse. For example, the contrasting compositions of the amphipathic lipid packing sensing (ALPS) motif of the ArfGAP1 protein (8) and the AH of α-synuclein (9) (Fig. 1, B and C) suggest that curvature sensing results from a delicate balance between the amino acid content on the apolar and polar sides of the helices (10,11). Moreover, it is important to note that some AHs (like α-synuclein) have a positive net charge, providing additional selectivity for anionic liposomes specifically (12). Together, the structural diversity among curvature sensors complicates reliable prediction of a given peptide's sensing ability simply from sequence-based physicochemical descriptors, like mean hydrophobicity 〈H〉, hydrophobic moment μ H (13), and net charge z.
Molecular dynamics (MD) simulations are a valuable asset in expanding our understanding of curvature sensing, since they can access the necessary molecular resolution that many experimental methods lack (14)(15)(16). To reduce system size and, consequently, reduce the computational cost, curved membranes are often represented as stretched flat membranes in MD simulations (Fig. 1A) (17,18), such that the lipid packing defects on the surface are similar and the consequent relative binding free energies correlate (19). This approximation is valid since biologically relevant bilayer curvatures are negligible on the length scale of peptides (≈5 nm).
Recently, we developed a highly efficient method that can quantify the relative free energy (ΔΔF) of curvature sensing by surface peptides (19). By redefining ΔΔF as a mechanical property, we realized that differential binding can be reinterpreted as the reduction of work required to stretch a membrane leaflet when the peptide is bound to it. In other words, a peptide that senses leaflet tension/curvature will also tend to induce leaflet tension/curvature, and therefore, those properties are two sides of the same coin. In line with this notion, many AHs have been shown to bind to flat membranes and then actively generate positive membrane curvature [e.g., the AHs of Epsin (20) and the N-BAR domain (21)].
A still unresolved key question is what physical characteristics differentiate AHs that specifically bind to curved membranes ("sensors") from AHs that bind without curvature specificity ("binders"). To date, this question has mostly been addressed by making strategically chosen point mutations in specific example cases (9,11,(22)(23)(24), but a fundamental thermodynamic understanding on how to distinguish sensors from binders among multiple chemically diverse classes of AHs is, to the best of our knowledge, still lacking. Here, we combine an evolutionary algorithm with coarse-grained MD simulations (Evo-MD) to design α-helical peptides that optimally recognize the curvature of lipid membranes, completely from scratch. Our main goal is to demonstrate how the unique synergy between Evo-MD and neural networks can enable the identification of the two major curvature sensing protein families (ALPS and α-synuclein) and their known mutants without using any of the available experimental input data, i.e., our approach is purely physics-based. Furthermore, we will illustrate that sensing and binding are phenomena that lie on the same thermodynamic continuum, with only subtle but explainable differences in membrane binding free energy, consistent with the serendipitous discovery of curvature sensors.

Designing the optimal curvature sensor
Evo-MD is a physics-based inverse design method that embeds MD simulations in a genetic algorithm (GA) framework (25). GAs are inspired by Darwinian evolution and can serve as a powerful tool for optimization problems in large discrete search spaces, like the 20 L possible peptide sequences (for 20 natural amino acids and peptide length L). Starting from a, in our case, random, initial subset ("population"), a GA iteratively (i) evaluates the desired property ("fitness") of the candidate solutions in the population, (ii) selects the best candidates as the "parents" for the next generation, and then (iii) performs genetic operations, like cross-over recombination and random point mutations to (iv) generate the next population ( Fig. 2A). While evolution proceeds, the population's average fitness will increase until it converges to an optimum. To date, GAs have mainly been applied to peptide optimization problems that involve protein-peptide interactions and use fitness functions based on physicochemical descriptors or information from databases (26)(27)(28)(29). In contrast, the fitness calculation in Evo-MD is based on ensemble averaging from (coarse-grained) MD simulation trajectories and is therefore completely data independent. In this physics-based approach, experimental data contribute to solving the optimization problem via the parametrization of the force field that is used in the simulations, the Martini model (30) in this study. Therefore, the main advantage is that Evo-MD will generate curvature sensing peptides without requiring any knowledge of existing curvature sensing peptides, of which too few examples exist to properly train a data-informed model. In addition, as opposed to data-trained models, physics-based inverse design does not tend to generate molecules that are (too) similar to the input data (31). In contrast, Evo-MD will search for a predefined thermodynamic optimum of sensing and generate physically optimal sequences that actually may differ from the biological optimum due to additional evolutionary constraints imposed by nature's complexity (e.g., solubility, protein-protein interactions, and trafficking).
The direction of simulated evolution by GAs is governed by the definition of the fitness function (the desired property). For the optimization of curvature sensing peptides, we aim to maximize the magnitude of the curvature sensing free energy ΔΔF that we can efficiently quantify using aforementioned mechanical free energy method (19). Our fitness function is the product of the ΔΔF value and a scaling factor c that equals 1 when located on the membrane surface and goes to 0 for transmembrane or soluble configurations (section S1). We emphasize that ΔΔF characterizes the relative affinity for lipid packing defects or equivalently positive leaflet curvature, analogous to the curvature-dependent binding constants (free energy of partitioning) measured in experimental model liposome assays (11,(22)(23)(24)(32)(33)(34)(35)(36) and analogous to the measured differences in the concentration of peripheral membrane proteins due to curvature-driven sorting in micropipette aspiration assays (37,38).
Since membrane-surface peptides are, in most cases, α-helical and the Martini force-field is unable to model protein folding events, we assume and fix helical secondary structure when generating the starting conformations for our peptides. To this end, and to reduce the search space, we excluded 10 amino acids with low αhelical propensities (39) (P, G, D, R, C, T, N, H, V, and I) while ensuring that every chemical subtype is represented in our final subset (comprising A, L, M, K, Q, E, W, S, Y, and F). We chose a fixed peptide length of 24 residues, which is in the typical range for curvature sensing peptides. Consequently, our search space contains 10 24 peptide sequences.
In the three independent Evo-MD runs we performed, we observed convergence within 25 iterations with the best candidates having a ΔΔF of around −32 kJ mol −1 (Fig. 2B). The consensus sequence logos (40) of the final generations show a strong enrichment of bulky hydrophobic residues, mainly F and W (Fig. 2, C to E, and movies S1 to S3). We can understand this result by returning to our earlier statement that "a peptide that senses tension/curvature will also tend to induce tension/curvature." Such induction of tension and concomitant membrane curvature in a membrane leaflet occurs via shallow insertions within the hydrophobic interior of the lipid membrane, i.e., the region directly below the head groups. We observe that the optimal sequences-the point of maximum leaflet tension generation with respect to the helix's central axis-is characterized by an insertion of 1.69 ± 0.05 nm from the bilayer center in a tensionless membrane (d mem in Fig. 2I) or, alternatively, 0.24 ± 0.05 nm below the average position of the phosphate groups (d PO4 in Fig. 2I). This is in quantitative agreement with predictions from membrane elastic theory, which suggest an optimal insertion of about 1.7 nm from the membrane center plane (41). Furthermore, the bulkier the peptide is, and thus the larger its excluded volume and effective helical radius, the more pronounced the induced leaflet tension will be.
Besides the abundant hydrophobic residues, we observed that the solutions of all three Evo-MD runs feature two charged residues (E or K; see Fig. 2, C to E). This numerical conservation of two charged amino acids suggests that this is the bare minimum of polar content that is necessary to maintain a surface orientation for such hydrophobic peptides (i.e., scaling factor c → 1). The sign of the charge appears to be irrelevant for the zwitterionic phosphatidylcholine (POPC) membranes we used here. Also, the exact position of these residues seems arbitrary, as long as the two charged residues end up on the same side of the helix in the folded conformation (Fig. 2, F to H).
The fact that all three randomly initiated Evo-MD runs produced peptide sequences with identical physical characteristics within the same number of iterations strongly suggests that this is indeed the global and not a local optimum. To probe the effect of only using the 10 most helix-prone amino acids, we performed an additional Evo-MD run with all 20 natural residues included. This, again, yielded peptides with the same physical characteristics but showed slower convergence (40 iterations) and higher diversity due to the vastly increased search space ( fig. S3). What this simulated evolution shows is that the GA has successfully selected a key aspect in curvature sensing, insertion of hydrophobic residues (22,36), which is then maximally amplified and exploited until the fitness converges. To such extent even, that the optimal "sensor" is so hydrophobic that it would likely stick to any membrane, regardless of curvature, thus being classified as a "binder" instead. What is immediately clear is that our optimized peptides strongly differ from the naturally evolved optima (e.g., the ALPS motif and α-synuclein), both in terms of ΔΔF (Fig. 2B) and in their chemical compositions (e.g., compare Fig. 1, B and C, with Fig. 2, F to H). Thus, our physics-based inverse design indicates that the distinction between curvature sensors and membrane binders can be considered as a continuum with a soft, subtle threshold at a relative binding free energy that is much lower than the theoretical optimum.
Biologically, the differences between the simulated optimum and naturally evolved peptides can be explained by considering the many boundary conditions imposed by the complex environment of in vivo systems. One of the most obvious and fundamental requirements is that proteins should be soluble in physiological buffer. The extremely hydrophobic GA-generated optima clearly fail this criterion and will readily aggregate and precipitate out of the solution. Our finding that curvature sensing is a hydrophobically driven process indicates that natural curvature sensing helical motifs have likely evolved from a trade-off between maximizing hydrophobic insertion (with one face of the helix) while remaining generally water soluble (the other face of the helix), giving rise to their amphipathic character. Also, since curvature sensing implies tension generation, peptides with a high |ΔΔF| could harm the integrity or shape of the membranes they adhere to. To demonstrate this, we performed additional simulations that show that the hydrophobic Evo-MD optimum is sufficiently tension inducing to generate positive curvature in a flat POPC membrane ( fig. S4). To circumvent such disruptive effects, an evolutionary pressure to limit this potency must exist. A neural network model to predict curvature sensing As a valuable byproduct of the iterative optimization process by Evo-MD, we obtained a large database of ≈54,000 unique sequences (all 24 residues long) and their respective sensing free energies (ΔΔF ) as calculated by MD simulations. With this wealth of data, we set out to train a convolutional neural network (CNN) that is able to predict curvature sensing ability from peptide sequence information only.
To enable the model to handle peptides shorter than 24 amino acids as well, we split the sequences in the original dataset at a random position, such that the resulting two fragments were at least seven residues long. Next, since ΔΔF depends linearly on length ( fig. S5), we interpolated the ΔΔF values for the split sequences, hereby tripling the dataset to ≈138,000 sequences (after discarding duplicate fragments). We refrained from extrapolating to sequences longer than 24 residues, since this would require additional assumptions on amino acid composition and potentially involve more complex tertiary structures that are inaccurately modeled by the Martini force field. A detailed description of the final training data is included in section S6.
As described previously in the context of activity prediction of helical antimicrobial peptides (42), we used one-hot encoded and zero-padded representations for the input sequences. These are then fed to two consecutive convolutional layers with max pooling, followed by a fully connected layer and a single output neuron to translate the connection weights into a float value: the predicted ΔΔF ( Fig. 3A; see section S7 for details on the optimization of the architecture and hyperparameters).
During the CNN training, minimization of the mean square error (MSE) converged after 18 epochs (table S3) to a MSE of 1.84 kJ 2 mol −2 for the validation set (25% random sample from the full data). For this validation set, we achieved excellent correlation (R 2 = 0.97) between the predicted and MD-calculated ΔΔFs (Fig. 3B). However, because sequences from the late iterations of the same Evo-MD run can be highly similar, the validation and training sets are arguably not fully independent. Therefore, to ultimately test our model, we predicted ΔΔF for 988 randomly generated sequences (between 7 and 24 residues long) that were not part of the training data and obtained a MSE of 2.92 kJ 2 mol −2 and a R 2 value of 0.66 when comparing the predicted values to ΔΔF calculated by MD simulations (Fig. 3C).
The trained neural network and all datasets are accessible at github.com/nvanhilten/CNN_curvature_sensing. Please note that the model should only be used for sequences between 7 and 24 amino acids long and that it assumes α-helical folding (as we did in the training data). On the basis of the performance of our model on the randomly generated test set, the root MSE of its predictions is ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi MSE p ¼ 1:71 kJ mol −1 , which is comparable to the typical errors obtained when calculating ΔΔF by MD simulation (e.g., compare the error bars in Fig. 4B).

Distinguishing sensing from binding
Now, with the MD quantification and CNN tools in hand, we can return to the key question posed in Introduction, namely, which, if any, characteristics can help us distinguish curvature sensors from binders, and what can relative binding free energies, like ΔΔF, teach us in this regard? To address this question, we composed a benchmark set of natural curvature sensing peptides (table S4), also including mutated variants that were empirically categorized as "nonbinders" (i.e., no affinity for any membrane) or binders (i.e., binding to membranes without curvature specificity). We should acknowledge the expert help of B. Antonny and R. Gautier in composing this list. In the following section, we will show that these benchmark peptides can be reliably categorized into separate regimes (see Fig. 4A) within a thermodynamic continuum, that is, of a continuous scale in terms of membrane binding free energies (ΔF) but features a sharp switch-like transition in terms of membrane partitioning behavior.
To fairly compare the sequences, we propose two correction factors to obtain an adjusted relative binding free energy ΔΔF adj . First, we linearly extrapolate the ΔΔF values of shorter peptides to their corresponding free energies if they were 24 residues long (ΔΔF L=24 , see section S5). Second, we realized that many of the peptides are cationic to improve interaction with (curved) anionic membranes that are abundant in nature. Since our MD simulations were performed with neutral POPC membranes, we hypothesized that the relative binding free energies would, in these cases, be underestimated and thus require a correction term c z z to account for this To determine the magnitude of c z (the relative free energy contribution per unit charge z), we performed additional MD simulations with anionic membranes [75% POPC and 25% 16:0-18:1 phosphatidylglycerol (POPG)] and found elevated relative binding free energies (|ΔΔF L=24 |), especially for the cationic peptides (table  S5). From the average difference between the ΔΔFs calculated on the different membranes, we obtained c z = −0.93 ± 0.89 kJ mol −1 per unit charge.
We calculated ΔΔF adj for the benchmark peptides with the CNN model and MD simulations (Fig. 4B). When ranking the peptides accordingly, we find that we can roughly reproduce the empirical qualification of nonbinders at the lower end and binders at the higher end of the ranking, as well as sensors in the middle. We find that the CNN-predicted ranking is in better agreement with the experimental trends than the results from MD simulations. We speculate that this is likely due to a smoothening effect, i.e., predictions of the MD simulations for individual peptide sequences are fully independent, whereas the CNN introduces effective correlations between sequences. Because the CNN is trained on the ensemble averaged values from many thousands of independent MD trajectories for different peptide sequences, we argue that disturbances ("noise") in chemical space (point mutations) and in the MD itself (limited sampling) are therefore smoothened out to such extent that the experimental trends are more robustly reproduced. Hence, we use the CNN-predicted ΔΔF adj values for the remainder of this paper.

Thermodynamic model of the sensing → binding transition
Along the lines of our current definitions, every peptide undergoing hydrophobically driven membrane binding is able to sense positive membrane curvature due to the increase in surface hydrophobicity upon bending. In other words, positive curvature enhances hydrophobically driven membrane binding. Empirically, however, a peptide is only classified as a curvature sensor if it only binds to membranes characterized by a high positive curvature. In this section, we will argue that the empirical classification of sensors versus binders can be intuitively understood from the population statistics of a two-state partition function.
Here, we define the following two states: (i) state m: The peptide is bound to the membrane; (ii) state s: the peptide is in the solution. We define the partitioning free energy difference between the two states as ΔF sm : the free energy of membrane binding minus the free energy of solvation with respect to the peptide in the gas phase. Using thermodynamic integration, we calculated ΔF sm for a nonbinder, a sensor, a binder, and the extremely hydrophobic Evo-MD optimum (Fig. 2, D and G) binding to a flat tensionless membrane (R = ∞; see section S10) and found that it linearly relates to ΔΔF adj (Fig. 4C). The reason for this linear correlation is that both membrane binding (ΔF sm ) and curvature sensing (ΔΔF adj ) are driven by hydrophobic interactions. Since we are interested in peptides that bind to curved membranes (10 ≤ R ≤ 100), we generalize ΔF sm to a vesicle radius-dependent form ΔF sm (R) (full derivation in section S10) and take a typical radius of R = 50 nm for the following calculations. At thermal equilibrium, the relative probability for a peptide to bind to a membrane P m is Boltzmann distributed. Consequently, P m is given by Here, N s represents the number of realizations within the accessible solvent, and N m represents the number of membrane binding realizations. The mathematical form of Eq. 2 resembles that of a socalled Fermi-Dirac function. If N s N m ¼ 1, this function features a sharp but continuous transition at ΔF sm = 0 (dashed curve in Fig. 4E), i.e., the point where membrane binding and peptide solubility are precisely in balance (P m = 0.5). However, the number of realizations in solution is expected to be much larger than the number of realizations associated with membrane binding, N s N m � 1, at typical lipid concentrations. Consequently, the transition point shifts to the right, i.e., peptide-membrane binding is favored over peptide solvation at the transition point P m = 0.5. The steep nature of this "switch function" notably explains why empirically classified sensing behavior for one peptide switches to binding behavior for another, possibly highly similar, peptide based on only subtle differences in (relative) membrane binding free energy, like we observed in this work (Fig. 4B). To finalize our model, we estimate the prefactor N s N m . The accessible membrane area in 1 liter of solution is A ¼ 1 2 cN A A lip , with c being the lipid concentration, N A being the Avogadro constant, and A lip being the area per lipid. The characteristic surface area of a helical peptide is roughly A p = 5 × 1 nm 2 , and, equivalently, its volume in solution is V p = 5 × 1 × 1 nm 3 . Taking c = 1 mM [a typical value in the concentration range used in experiments (11,(22)(23)(24)(32)(33)(34)(35)(36)] and A lip = 0.64 nm 2 , we obtain N s N m ¼ V=V p A=A p ¼ 5:2 � 10 3 . When we plug this number into Eq. 2, we find that the free energy of membrane binding outperforms the free energy of solvation by −21.2 kJ mol −1 at the transition point P m = 0.5 (Fig. 4E), i.e., the transition point is associated with favorable but relatively weak membrane binding. Because of the sharp transition behavior, sensors are envisioned as peptides with a ΔF sm value near the transition point (0.05 ≤ P m ≤ 0.95; the orange area in Fig. 4E). According to our model, a peptide with P m < 0.05 would be classified as a "nonbinder," and a peptide with P m > 0.95 would be deemed a binder. When overlaying the values for the 19 benchmark peptides we introduced earlier (vertical lines in Fig. 4E and table S6) and considering a typical vesicle radius R = 50 nm, we find that 3 of 3 nonbinders, 8 of 12 sensors, and 3 of 4 binders are classified correctly, i.e., in agreement with the original experiments. Such categorization would have been impossible with crude physicochemical descriptors like mean hydrophobicity 〈H〉 or hydrophobic moment μ H despite them (weakly) correlating with ΔF sm (section S12).
We defined the sensing regime rather generously (0.05 ≤ P m ≤ 0.95) to accommodate for the fact that both the empirical classification and the positioning of the Fermi-Dirac curve are sensitive to the lipid concentration c, which can greatly vary between experiments. Moreover, imperfect helical folding and differing lipid compositions, cell systems, and read-out types can complicate the direct comparison between empirical classifications and our computational models. In terms of the curved membrane binding free energy ΔF sm (R = 50), the lower and upper boundaries of sensing (0.05 ≤ P m ≤ 0.95) correspond to −13.9 and −28.5 kJ mol −1 , respectively (orange area in Fig. 4E). In other words, when the work required to pull a peptide from the membrane of a vesicle with radius R = 50 is between those values, it can be classified as a sensor. In terms of the length-and charge-adjusted relative binding free energy (ΔΔF adj ), i.e., the preference for highly curved membranes [vesicles with radius R = 12.5 nm (19)], sensors fall between −6.4 and −10.0 kJ mol −1 (orange area in Fig. 4D).
As an example, we highlighted four variants of ALPS GMAP-210 (23,24) for which the relative differences in experimental sensing/ binding behavior were correctly captured (Fig. 4D). This is a notable example because the sequences are similar: The only difference between the original sensor ALPS GMAP-210 and the nonbinding variant ALPS GMAP-210 (L12D) is a single point mutation (L → D) that disturbs the peptide's hydrophobic face. As can be expected, and in line with the experimental findings, the inverse sequence ALPS GMAP-210 (inv) scores the same as the original peptide and is thus categorized as a sensor as well. Last, the condensed version ALPS GMAP-210 (cond) was correctly identified as a binder. With this highlighted example, we demonstrate that our method and thermodynamic model can (i) pick out features as subtle as single point mutations, (ii) resolve the resulting differences in relative free energy, and (iii) correctly categorize the consequent sensing behavior.
We also included the antiviral peptide hepatitis C virus (HCV) AH, which specifically ruptures vesicles with a high curvature (e.g., small enveloped viruses) (3, 4) (black dotted line in Fig. 4E). To date, HCV AH is the only example of a clinically relevant curvature selective antiviral peptide (5). When we plug in the previously calculated free energy value (19) for HCV AH, we find that this peptide falls into the binder regime. This is consistent with evidence that the vesicle size specificity of this peptide is due to curvature-specific pore formation and not to curvature-specific binding (i.e., curvature sensing) (43). After all, subtle binding may not be optimal for pore formation, since the subsequent induction of tension should be sufficient to rupture the membrane.
Following the evaluation of HCV AH in this regard, we argue that the most promising range to find potent curvature-specific antiviral agents is therefore near the transition zone between sensing and binding (i.e., P m → 1), since these peptides (i) may still benefit from some "curvature sensing" (predominant binding to higher curvatures) and (ii) pack a larger punch than biological sensors in terms of meeting the tension induction threshold necessary to deform/perforate viral membranes but are (iii) not so potent that they also rupture the host cell membrane. The latter is helped by the fact that the host membrane is likely more resilient to the disruptive actions of peptides than the viral membrane due to membrane-stabilizing proteins and active feedback mechanisms. This means that a considerable part of the selectivity of membrane-targeting drugs is likely due to a difference in drug resistance (membrane resilience) rather than actual differential binding.
Last, we plotted the membrane binding probability P m as a function of vesicle radius R for different ΔΔF values (Fig. 4F) and find that this is precisely in line with the cartoon in Fig. 4A. For nonbinders (ΔΔF > −6.4 kJ mol −1 ), membrane binding is rare (low P m ) for biologically relevant vesicle sizes (10 ≤ R ≤ 100). Conversely, binders (ΔΔF < −10.0 kJ mol −1 ) have a high membrane binding probability (P m → 1) regardless of the vesicle radius. It is exactly the sensor region in between that defines the area where peptides can display curvature selectivity, i.e., a wide range of probabilities for a wide range of vesicle radii.

DISCUSSION
We have illustrated the utility of a physics-based generative model (Evo-MD) to explore and simultaneously rationalize the mechanisms of how peptides sense membrane curvature. Initially, we set out to optimize curvature sensing by resolving the sequence that maximizes the relative affinity for lipid packing defects. Instead, we ended up with the optimal binder. This finding led to the important realization that curvature sensing and membrane binding are phenomena that lie on the same thermodynamic continuum (Fig. 4E).
Naturally evolved curvature sensors, such as the ALPS motif and α-synuclein, are chemically diverse but turn out to be remarkably similar in terms of partitioning free energies, which explains their functional similarities. In this work, we described the thermodynamic regime that defines the curvature sensing behavior of peptides. Given how narrow this energetic "sensing window" is, it is expected that the discovery and design of curvature sensors has, to date, been rather serendipitous. Having identified this "window of opportunity" in terms of relative binding free energy can facilitate the discovery of previously unidentified curvature sensing peptides since we now know where to look for them. The existence of the here-resolved thermodynamic sensing regime can also be intuitively understood from an evolutionary biological perspective. Curvature sensing motifs within naturally evolved proteins must fulfill the following two criteria: They should (i) predominantly bind to curved membranes and (ii) conserve the structural integrity of the membranes they adhere to. The here-observed linear correlation between sensing and overall membrane binding (Fig. 4C) dictates that these criteria are only met in the weak binding regime. These arguments are all in full agreement with the earlier hypothesis that curvature sensing in nature is a subtle balance between overall membrane binding and specific curvature recognition (2). Also, in this weak binding regime, the small leaflet strain induced by peptides is able to facilitate a largely inert and thus biologically functional sensing phenotype that does not easily lead to membrane rupture/deformation.
We demonstrated a fruitful synergy between a physics-based generative model (Evo-MD) and a CNN, which not only markedly accelerated the high-throughput evaluation of peptide sequences but also improved the accuracy of prediction compared to the original molecular simulations. It is important to stress the key role of Evo-MD, in that it yields sequences over the whole range of ΔΔF by gradually maximizing the relevant chemical property in a wellspaced manner, in our example, even up to the thermodynamic optimum. Using these data to train a neural network model has the important advantage that it encompasses the full thermodynamic range of possibilities over a vast search space of 20 24 sequences, whereas a dataset of natural peptides (if available in the first place) would be strongly constrained to a certain biologically feasible regime that only comprises peptides with highly similar physicochemical characteristics. We argue that training data generated with Evo-MD can therefore substantially improve both the applicability domain and the accuracy of neural network models, despite many of the generated sequences not being necessarily biologically relevant. This principle is equivalent to fitting an unknown function to data points that are well spaced over the whole range of the applicability domain versus data points that are only clustered within a narrow window. In particular, precise knowledge of the maxima (and minima) of a function-which a physics-based optimization resolves-will benefit the quality of a fit or model, also within the biologically relevant domain of the search space. We postulate that a subsequent restriction of the search space within or near the hereresolved sensing regime can enhance the discovery of curvature sensing motifs in natural proteins, as well as their de novo generation.
Last, we envision an important potential application in the computational design of peptide sensors that recognize membranes with other aberrant characteristics, such as a distinct lipid composition (e.g., bacteria and cancer cells). Since (selective) membrane binding results in the generation of leaflet tension, membrane binding peptides have an inherent membrane-destabilizing propensity and the ability to lower the energetic cost of the highly curved interface of toroidal pores. This is particularly the case for the hydrophobically driven membrane binding peptides we discussed in this work. The simultaneous encoding of selective membrane binding and an active drug mode, such as the induction of membrane lysis, is therefore a realistic avenue to explore further. An important advantage of physics-based generative models over existing data sciencebased generative models here is their unique ability to systematically explore the drug therapeutic potential of distinct relative binding (ΔΔF ) regimes by restricting the generation of peptide sequences within predefined boundary values of ΔΔF, for example, via the straightforward introduction of a bias/constraint to the fitness function. This can enable the targeted exploration of different "windows of opportunity" similar to aiming a gun at different targets.

Evo-MD
The Evo-MD code was adapted from previous work (25) to include the preparation, production, and fitness evaluation modules that are specific to the curvature sensing problem. The code is available at github.com/nvanhilten/Evo-MD_curvature_sensing.
Every candidate solution is evaluated by a fitness function that is based on a free energy calculation from coarse-grained MD simulations, using the Martini 3 force field (30). The relative binding free energy (ΔΔF ) is calculated from two POPC membrane simulations at a constant area: one at tensionless conditions [σ(A 0 )] and one at stretched conditions [σ(A * )] with a relative strain of ε ¼ A � À A 0 A 0 ¼ 0:165, such that the outer leaflet bares lipid packing defects similar to the outside of a highly curved liposome (radius R = 12.5 nm). This method is described in detail in previous work (19) and in section S1.
To focus the optimization toward surface peptides, we applied a cosine scaling factor c (Eq. 3 and section S1) that is 1 when the peptide is at the membrane surface (d mem = b) and 0 when the peptide is in solution (d mem ≥ 2b) or fully adopts an intermembrane or transmembrane configuration (d mem ≤ 0). The reference values for b are based on the monolayer thickness of a tensionless POPC membrane [b = 1.90 nm for σ(A 0 )] and a stretched POPC membrane [b = 1.79 nm for σ(A * )] (19) within the Martini 3 force field. The lowest value for c (tensionless or stretched condition) is used in the fitness scaling Evo-MD production runs were performed with a population size of 144 peptides, from which a parent pool of the "fittest" 36 sequences was selected. Details on tuning these settings are described in section S2. From the parent pool, two parents were selected randomly and combined through one-point crossover recombination: The two parent sequences are split at the same random position i and the tail ends (i → 24 − i) are swapped. Next, for every position in the resulting two "children" peptides, random point mutation is applied with the probability P mut ¼ 1 L , with sequence length L = 24, such that every sequence bares, on average, one point mutation.
To retain highly scoring peptides, the best candidate solution of every generation ("fitness elite") was copied to the next. In addition, sequences that were rerun more than three times already were also retained ("rerun elites"). The fitness values for these elite sequences were updated after every rerun with the average value of the fitness values from all individual iterations. In this way, highly scoring peptides get increasingly better sampling throughout the evolution. Simulation details All MD simulations were performed with GROMACS 2019.3 (44) with fixed x-and y-dimensions (constant area), such that the "low tension" membrane was tensionless (ϵ = 0, A = A 0 ), and the "high tension" was at a 16.5% relative strain (ϵ = 0.165, A = A * ) (19). Berendsen pressure coupling (45) was applied only in the z-dimension (1-bar reference pressure with a 4.5 × 10 −5 bar −1 compressibility, τ P = 6 ps). A constant temperature of 310 K was maintained by the velocity rescaling thermostat (46) (τ T = 1 ps). We used the coarse-grained Martini force field, version 3.0.0 (30), at a 30-fs time step. Van der Waals interactions were calculated with the shifted Verlet cutoff scheme (47), and reaction-field electrostatics (48) describe the coulomb potentials, both with a 1.1-nm cutoff. The neighbor list was updated every 20 steps.
Within Evo-MD, the system setup was automated as described in previous work (19). Systems comprise 128 Martini 16:0 to 18:1 POPC lipids and 1800 Martini water beads. Atomistic helical peptide models were generated using PeptideBuilder (49) and then coarse-grained by martinize2, VerMoUTH (50). They were then inserted into the two preequilibrated membrane systems at d mem = 1.5 nm from the membrane center plane. For charged peptides, systems were neutralized by adding Na + and Cl − counter ions. Steepest descent minimization with soft-core potentials (0.75 coupling for the Van der Waals interactions) was performed to solve clashes. Within Evo-MD, a total run time of 400 ns was used for each simulation. All non-Evo-MD simulations (Figs. 3C and 4B, fig. S5, and table S5) were run for at least 800 ns.

Supplementary Materials
This PDF file includes: Sections S1 to S13 Figs. S1 to S8 Tables S1 to S6 Legends for movies S1 to S3 References Other Supplementary Material for this manuscript includes the following: Movies S1 to S3 View/request a protocol for this paper from Bio-protocol.