Peptide Vaccine Design by Evolutionary Multi-Objective Optimization

1 Introduction

Peptide vaccines, composed of a set of peptides with the ability to selectively activate and proliferate T cells, are becoming increasingly essential in clinical treatment of a variety of diseases, including human immunodeficiency virus [Kong et al., 2016; Xu et al., 2018], Alzheimer's disease [Davtyan et al., 2013; Wang et al., 2017], and various forms of cancer [Malonis et al., 2020; Abd-Aziz and Poh, 2022]. Compared to traditional live-attenuated vaccines, peptide vaccine design offers a high degree of specificity and avoids the risks of infection caused by using entire pathogens.

Epitopes on peptides can bind to Major Histocompatibility Complex (MHC) molecules within the human body, thereby eliciting specific immune responses. Machine learning has greatly aided peptide (epitope) identification, involving detecting peptides that trigger immune responses [Jurtz et al., 2017; Reynisson et al., 2020]. Many impressive works have benefited from advanced peptide prediction methods [Baruah and Bose, 2020; Bhattacharya et al., 2020]. After predicting pathogen peptides, the main challenge in peptide vaccine design is selecting an effective peptide subset. This task is complicated by the extensive allelic diversity of human MHC molecules, resulting in significant variability in peptideMHC combinations across individuals [Robinson et al., 2020; Barker et al., 2023]. Therefore, selecting a diverse subset of peptides that can be universally effective for a broad population poses a considerable challenge.

1.1 RelatedWork

Most previous work on peptide vaccine design typically relies on discrete optimization methods, such as integer linear programming and genetic algorithms, to maximize population coverage [Oyarzun and Kobe, 2015]. Liu et al. [2022] emphasized the importance of incorporating redundancy to boost the likelihood of effective immunogenic peptide presentation. They treated peptide vaccine design as a maximum N-times coverage problem, aiming to select peptides that ensure at least N immunogenic peptide-MHC interactions per individual, and solved it with heuristic methods. However, this problem is NP-hard and, being non-submodular, does not allow polynomial-time constant factor approximation. This formulation also tends to underestimate the utility of a peptide until N-times coverage is achieved.

Dai et al. [2023] formulated peptide vaccine design as a constrained optimization problem, aiming to maximize the expected number of peptide-MHC bindings across a broad range of populations while choosing a limited set of diverse peptides. They introduced constraints on the size of the peptide subset and promoted dissimilar redundancies to avoid correlated failures, leading to the cardinality and pairwise constraints. The objective function, i.e., the expected number of peptide-MHC bindings across a wide population range, is proved to be both monotone and submodular, which enables their greedy algorithm, Optivax-P, to provide an approximation guarantee of max{OPT/(1 + ∆), OPT2/2}, where ∆ > 0, OPT denotes the optimal function value, and OPT2 is the optimal function value under additional constraints.

1.2 Our Contribution

Considering that the performance of the greedy algorithm may be limited due to its greedy nature, this paper proposes a new framework based on Evolutionary Multi-objective Optimization [Knowles et al., 2001; Friedrich and Neumann, 2015; Qian et al., 2015] for Peptide Vaccine Design, briefly called PVD-EMO, which reformulates peptide vaccine design as a bi-objective optimization problem that maximizes the expected number of peptide-MHC bindings and minimizes the number of selected peptides simultaneously, and employs a Multi-Objective Evolutionary Algorithm (MOEA) to solve it. PVD-EMO can be equipped with any MOEA to solve this bi-objective problem, and we employ the theoretically grounded GSEMO [Laumanns et al., 2004] as well as the popular NSGA-II [Deb et al., 2002]. To boost efficiency, we also incorporate warm-start and repair strategies into MOEAs. The warm-start strategy seeds an initial population of PVD-EMO with a range of solution sizes, including the output solution of the previous greedy algorithm OptivaxP [Dai and Gifford, 2023]. Meanwhile, the repair strategy is designed to steer PVD-EMO towards feasible solutions when it encounters infeasible regions, thereby improving its exploratory potential. We prove that the warm-start strategy ensures that PVD-EMO maintains the same worst-case approximation guarantee as the previous greedy algorithm OptivaxP. Additionally, by using an illustrative example of peptide vaccine design, we prove that Optivax-P will get trapped in local optima, while PVD-EMO can help avoid them. Experiments on a peptide vaccine design for COVID-19, caused by the SARS-CoV-2 virus, demonstrate the superiority of PVDEMO over the state-of-the art algorithm Optivax-P.

2 Peptide Vaccine Design

Let V = {v1, v2, . . . , vn} denote the set of peptides, and M denote the set of MHC genotypes observed in the population. Peptide vaccines work by selecting an effective subset of peptides displayed on the cell surface of MHC proteins. The resulting peptide-MHC complexes activate the cellular immune system. Let pv,m be the probability that a peptide v ∈ V is displayed by an individual's MHC genotype m ∈ M, with the assumption that these probabilities are independent across different peptides. The occurrence of peptide v being displayed is termed a peptide-MHC hit, represented by the indicator function I(pv,m), which equals 1 iff the event occurs.

Different MHC alleles have different peptide binding properties, so it is important for a vaccine to trigger multiple peptide-MHC bindings, ensuring redundancy in the activation of T cell clonotypes in an individual. Given a subset S of selected peptides, for an MHC genotype m ∈ M, the number of peptides in S that are displayed can be represented as ∑v∈S I(pv,m). Since too much redundancy may lead to unnecessary burden on the immune system, a threshold parameter N ≥ 0 is used to limit the additional benefits until a person attains the desired N hits, i.e., min{∑v∈S I(pv,m),N}.

Thus, the expected number of peptide-MHC bindings for a subset S across the whole MHC genotypes set M is

where w(m) denotes the weight corresponding to the percentage of the population with the genotype m.

To obtain f(S), we calculate E[min{∑v∈S I(pv,m),N}]

for each m ∈ M, and then sum up these expectations across all m. Assume that the indicator variables I(pv,m) are independent for each MHC genotype m ∈ M. For a given MHC genotype m, let Si ⊆ S denote the subset consisting of the first i peptides of S (assume an arbitrary order of the peptides in S), and Yi = ∑v∈Si I(pv,m) denote the sum of |Si| independent Bernoulli trials. The distribution of Y0 is trivially P(Y0 = 0) = 1, and we can iteratively compute the distribution of Yi as P (Yi = y) = P (Yi-1 = y - 1) · p + P (Yi-1 = y) · (1-p), where p denotes the probability pv,m that the new peptide v ∈ Si \ Si-1 is displayed by the MHC genotype m. Then, the distribution of Z = min{∑v∈S I(pv,m), N} satisfies that ∀z < N : P(Z = z) = P(Y|S| = z), and P(Z = N) = P(Y|S| ≥ N), where Y|S| is just ∑v∈S I(pv,m). The expectation is then computed using the resulting distribution.

The objective function, given by Eq. (1), is monotone and submodular [Dai and Gifford, 2023]. Let R denote the set of reals. A set function f : 2V → R is monotone if ∀X ⊆ Y ⊆ V , f(X) ≤ f(Y ). As more peptides will not worsen the value, the monotonicity is satisfied naturally. A set function f is submodular [Nemhauser et al., 1978] if it satisfies the "diminishing returns" property, i.e., ∀X ⊆ Y ⊆ V, v /∈ Y ,

which implies that the benefit of adding a peptide to a set will not increase as the set extends.

The peptide vaccine design problem is subject to two types of constraints. The first is a cardinality constraint that ensures the selected peptide subset S not exceed a given size k, i.e., |S| ≤ k, which is crucial for increasing the vaccine's stability, as well as reducing production costs. The second type involves pairwise constraints, which prevent any two peptides in the subset S from being similar. Maintaining dissimilarity among peptides is crucial for vaccine design, as similar peptides might fail for the same reason, reducing the vaccine's overall efficacy. Let G = (V,E) denote a graph, where the vertices correspond to the peptide set V , and the edges connect peptides that are deemed similar. A solution that satisfies the set of pairwise constraints corresponds to an independent set in the graph G = (V,E); that is, for any two peptides vi, vj in the subset S ⊆ V , there is no edge (vi, vj) in E.

In [Dai and Gifford, 2023], peptide vaccine design has been formulated as the problem of maximizing the objective function f(S), which characterizes the expected number of peptide-MHC bindings across a broad range of populations, subject to a cardinality constraint and a set of pairwise constraints, as presented in Definition 1.

Definition 1 (Peptide Vaccine Design). Given a set of peptides V = {v1, v2, · · · , vn}, a set of MHC genotypes M, probabilities pv,m of binding between each peptide v and each MHC genotype m, a weight function w(m) corresponding to the percentage of the population with each MHC genotype m, a threshold N ≥ 0, a budget k, and a similarity graph G = (V,E), the goal of peptide vaccine design is to find a subset S ⊆ V of peptides that maximizes f(S) subject to a size constraint and a set of pairwise constraints. That is:

2.1 Previous Algorithm

For the peptide vaccine design problem in Definition 1, Dai and Gifford [2023] introduced a greedy algorithm, Optivax-P, which starts from an empty set and iteratively adds peptides with the highest marginal gain on f, satisfying given constraints. Let ∆ denote the maximum degree of the similarity graph G. When ∆ = 0, implying no pairwise constraints, the problem simplifies to optimizing a monotone submodular function with a cardinality constraint, where Optivax-P achieves the optimal (1 - 1/e) polynomial-time approximation ratio [Nemhauser et al., 1978]. For ∆ > 0, indicating pairwise constraints, Optivax-P can find a solution ˆ S satisfying f( ˆ S) ≥ max{OPT/(1 + ∆), OPT2/2}, where OPT and OPT2 denote the optimal function values for the original and a more constrained problem, respectively. For the more constrained problem, G is replaced by G2 , linking vertices within two steps in G.

3 PVD-EMO Framework

Inspired by the excellent performance of MOEAs for solving general subset selection problems [Friedrich and Neumann, 2015; Qian et al., 2015; Qian et al., 2017a; Qian et al., 2017b; Qian et al., 2018; Qian et al., 2019; Qian et al., 2020; Qian, 2020; Bian et al., 2021; Qian, 2021; Bian et al., 2022; Qian et al., 2022; Roostapour et al., 2022; Qian et al., 2023; Zhang et al., 2023], we propose a new Peptide Vaccine Design framework based on Evolutionary Multi-objective Optimization, called PVD-EMO. A subset S of V can be naturally represented by a Boolean vector s ∈ {0, 1}n , where the i-th bit si = 1 iff the i-th peptide in V is contained by S. We will not distinguish s ∈ {0, 1}n and its corresponding subset {vi ∈ V | si = 1} for notational convenience.

As presented in Algorithm 1, PVD-EMO first reformulates the original peptide vaccine design problem in Definition 1 as a bi-objective maximization problem

That is, the first objective f1 equals the original objective f (i.e., the expected number of peptide-MHC bindings across a broad range of populations) for feasible solutions satisfying the cardinality and pairwise constraints, while -1 for infeasible ones; the second objective f2(s) = -|s| = - ∑ni=1 si is the opposite of the subset size. The domination relationship in Definition 2 is used to compare solutions. A solution

Algorithm 1 PVD-EMO Framework

Input: a peptide vaccine design problem instance, a budget k and a similarity graph G = (V, E)

Output: a subset of V with size at most k

Process:

• 1: Construct two objective functions f1(s) and f2(s) to be maximized, as presented in Eq. (3);
• 2: Apply an MOEA to solve the bi-objective problem;
• 3: return the best feasible solution in the final population generated by the MOEA

is Pareto optimal if no other solution dominates it. The collection of objective vectors of all Pareto optimal solutions is called the Pareto front.

Definition 2 (Domination). For two solutions s and s′ ,

1. s weakly dominates s′ (i.e., s is better than s′ , denoted by s ⪰ s′ ) if ∀i : fi(s) ≥ fi(s′ );
2. s dominates s′ (i.e., s is strictly better than s′ , denoted by s ≻ s′ ) if s ⪰ s′ ∧ ∃i : fi(s) > fi(s′ );
3. s and s′ are incomparable if neither s ⪰ s′ nor s′ ⪰ s.

After constructing the bi-objective problem in Eq. (3), PVD-EMO employs an MOEA to solve it, as shown in line 2 of Algorithm 1. Evolutionary algorithms (EAs), inspired by Darwin's theory of evolution, are general-purpose randomized heuristic optimization algorithms [B¨ack, 1996], mimicking variational reproduction and natural selection. Starting from an initial population of solutions, EAs iteratively reproduce offspring solutions by crossover and mutation, and select better ones from the parent and offspring solutions to form the next population. The population-based search of EAs matches the requirement of multi-objective optimization, i.e., EAs can generate a set of Pareto optimal solutions by running only once. Thus, EAs have become the most popular tool for multi-objective optimization [Coello et al., 2007; Hong et al., 2021; Yang et al., 2024; Liang et al., 2024], and the corresponding algorithms are also called MOEAs. After running a number of iterations, the best feasible solution (i.e., the solution having the largest f value while satisfying the constraints) will be selected from the final population as the output, as shown in line 3 of Algorithm 1. Note that the aim of PVD-EMO is to find a good solution of the original peptide vaccine design problem in Definition 1, rather than the Pareto front of the reformulated bi-objective problem in Eq. (3). That is, the bi-objective reformulation is an intermediate process. The introduction of the second objective f2 can naturally bring a diverse population, which may lead to better optimization performance.

PVD-EMO can be equipped with any MOEA, and we employ the theoretically grounded GSEMO [Laumanns et al., 2004] as well as the popular NSGA-II [Deb et al., 2002]. To boost efficiency, we incorporate warm-start and repair strategies into MOEAs, and also design a strategy to accelerate the objective evaluation. Next, we will introduce them in detail.

3.1 PVD-GSEMO-WRAlgorithm

PVD-EMO which employs the GSEMO algorithm and uses both Warm-start and Repair strategies is called PVD-

Algorithm 2 PVD-GSEMO-WR Algorithm

Input: a peptide vaccine design problem instance, a budget k and a similarity graph G = (V, E)

• Output: a subset of V with size at most k Process: 1: Construct two objective functions f1(s) and f2(s) to be maximized, as presented in Eq. (3); 2: Initialize the population P byWarm-Start Strategy; 3: repeat 4: Choose s from P uniformly at random; 5: Create s′ by flipping each bit of s with prob. 1/n; 6: s′′ ←Repair Strategy (s, s′ , G); 7: if ∄z ∈ P such that z ≻ s′′ then 8: P ← (P \ {z ∈ P | s′′ ⪰ z}) ∪ {s′′ } 9: end if 10: until some criterion is met 11: return the best feasible solution in P

GSEMO-WR, as presented in Algorithm 2. PVD-GSEMOWR first constructs the bi-objective problem in Eq. (3) in line 1 of Algorithm 2. After that, it starts with an initial population created by the warm-start strategy in line 2, and iteratively improves the quality of solutions in the population P (lines 3-10). The warm-start strategy presented in Algorithm 3 randomly generates a feasible solution for each size in {0, 1, . . . , k - 1}, while for the solution of size k, it uses the output from the previous greedy algorithm, OptivaxP [Dai and Gifford, 2023]. Among these solutions, the nondominated ones will be included in the initial population P. In each iteration, a parent solution s is selected from the current population P uniformly at random (line 4), and used to generate an offspring solution s′ by bit-wise mutation (line 5), which flips each bit of s independently with probability 1/n. Then the offspring solution s′ will go through a repair strategy (line 6), as presented in Algorithm 4, which can fix an offspring solution s′ such that it no longer violates the pairwise constraints, and will be introduced later. The repaired offspring solution s′′ equals to s′ if s′ has already been feasible. Then, s′′ is used to update the population P (lines 7-9). If s′′ is not dominated by any solution in P (line 7), it will be added into P , and meanwhile, those solutions weakly dominated by s′′ will be deleted (line 8). This updating procedure makes the population P always contain incomparable solutions. Furthermore, P always contains feasible solutions, because 1) the repair strategy can fix the violation on the pairwise constraints; 2) a solution violating the size constraint has bad values on both objectives according to Eq. (3), and will not be included into the population.

In Algorithm 4, the repair strategy reviews each bit of parent solution p and offspring solution o (lines 1-7). When a bit i flips from 0 in p to 1 in o (line 2), it identifies a set Q of indices in o that are connected to peptide vi in edge set E (line 3). The algorithm then randomly keeps one index q from Q unchanged (line 4) and sets the other connected bits in Q\q to 0 (line 5). This ensures that no pairwise constraints are violated in the first i bits of o. After processing all bits, the repaired o, free from pairwise constraint violations, is returned (line 8).

Algorithm 3 Warm-Start Strategy

Input: solution sg output by the greedy algorithm Optivax-P Process:

1: P = {sg}; 2: for i = 0 to k -1 do 3: Randomly create a feasible solution s of size i; 4: if ∄z ∈ P such that z ≻ s then 5: P ← (P \ {z ∈ P | s ⪰ z}) ∪ {s} 6: end if 7: end for

• 8: return an initial population P

Algorithm 4 Repair Strategy

Input: a similarity graph G = (V,E), a parent solution p and an offspring solution o

Process:

• 1: for i = 1 to n do 2: if oi = 1 and pi = 0 then 3: Q = {j | (vi, vj) ∈ E & oj = 1} ∪ {i}; 4: Choose q fromQ uniformly at random; 5: Set op to 0, for any p ∈ Q \ {q} 6: end if 7: end for 8: return the offspring solution o

3.2 PVD-NSGA-II-WRAlgorithm

NSGA-II is a popular MOEA which incorporates two substantial features, i.e., non-dominated sorting and crowding distance. For a detailed description of NSGA-II, please refer to [Deb et al., 2002]. PVD-EMO which employs NSGA-II for multi-objective optimization and incorporates both Warmstart and Repair strategies is called PVD-NSGA-II-WR.

In Section 4, we will prove that PVD-GSEMO-WR and PVD-NSGA-II-WR achieve the same theoretical guarantee as Optivax-P [Dai and Gifford, 2023] while better avoiding local optima. In Section 5, we will show their superior performance over Optivax-P in real-world experiments.

3.3 Acceleration of Objective Evaluation

In each iteration of PVD-EMO, we need to evaluate the objective value of a newly generated solution s′′ , i.e., f(S′′ ) = ∑m∈M w(m)·E[min{∑v∈S′′ I(pv,m),N}], where S′′ is the corresponding subset of s′′ . For each m ∈ M, a random variable Y = ∑v∈S′′ I(pv,m) is the sumof |S′′ | ≤ k independent Bernoulli trials I(pv,m), and we can get the probability distribution of Y by iterated convolutions, which costs O(k2 ) time. E[min{∑v∈S′′ I(pv,m),N}] then can be calculated by using E[h(y)] = ∑y h(y) · P(Y = y), where h(y) = min{y, N}, which requires O(k3 ) time. Thus, the total time of computing f(S′′ ) is O(|M|k3 ), which is expensive.

Next, we give an acceleration strategy of evaluating f (S′′ ). For each m ∈ M, let DX denote the probability distribution of a random variable ∑v∈X I(pv,m), and DX j denote the probability of attaining j hits from a peptide subset X, i.e., DX j = P(DX = j). The probability distribution D∅

of an empty subset satisfies that D∅ 0 = 1 and D∅ j = 0 for all 0 < j ≤ k. To accelerate the evaluation process, we build a recursive relation between the distributions DS′′ and DS , where S is a parent solution of S′′ . Specifically, we can compute DS′′ based on DS by recursively adding peptides from X+ = {vi | vi /∈ S ∧ vi ∈ S′′ } and then recursively deleting peptides from X- = {vi | vi ∈ S ∧ vi /∈ S′′ }. The distribution DX∪{vi} is updated based on DX :

Similarly, the probability distribution DX\{vi} is updated by backtracking from DX :

Note that the commonly used bit-wise mutation operator in line 5 of Algorithm 2 will flip one bit in expectation, that is, it will add or delete one peptide in expectation. Consequently, the expected time complexity for calculating the probability distribution DS′′ of an offspring solution S′′ is O(k), which is faster than the O(k2 ) time required by direct computation. We store the utility function h(y) = min{y,N} in a data structure that supports random access, such as an array, which allows us to compute its expectation in O(k) time. By implementing these optimizations, the overall time required to compute the objective function can be reduced from O(|M|k3 ) to O(|M|k).

4 Theoretical Analysis

For peptide vaccine design in Definition 1, the objective function f(s) has been proved to be monotone and submodular [Dai and Gifford, 2023]. The greedy algorithm Optivax-P can provide an approximation guarantee of max{OPT/(1 + ∆), OPT2/2}, where ∆ > 0, OPT denotes the optimal function value of Eq. (2), and OPT2 is the optimal function value under additional constraints. By maximizing f1(s) = f(s) and f2(s) = -|s| simultaneously, PVD-GSEMO-WR or PVD-NSGA-II-WR can achieve the same theoretical guarantee as Optivax-P, as shown in Theorem 1.

Theorem 1. For peptide vaccine design in Definition 1, PVDGSEMO-WR, or PVD-NSGA-II-WR with a population size of at least 4(k + 1), can achieve the same approximation guarantee as the previous greedy algorithm Optivax-P.

The proof is straightforward, because PVD-GSEMO-WR and PVD-NSGA-II-WR begin with an initial population containing the greedy solution output by Optivax-P, due to the warm-start process. The population update mechanism of PVD-GSEMO-WR ensures that if a solution is deleted, there must be another solution in the population that weakly dominates the deleted one. According to the recent theoretical work [Zheng et al., 2022], if the population size is at least four times of the size of the Pareto front (which is no larger than k + 1 here), PVD-NSGA-II-WR will keep at least one solution for each objective vector in the first non-dominated front of the population. Consequently, PVD-GSEMO-WR and PVD-NSGA-II-WR are guaranteed to obtain a solution that is at least as good as the one obtained by Optivax-P.

Using a peptide vaccine design example from Definition 3, Theorem 2 shows that Optivax-P falls into local optima, while PVD-GSEMO-WR can avoid them and achieve the global optimum, which excludes the peptide v1. The proof shows that due to its greedy nature, Optivax-P first selects v1 and is misled. In contrast, PVD-GSEMO-WR either avoids v1 via bitwise mutation or escapes local optima using a repair strategy, finally following the population's path to the global optimum.

Definition 3 (An Example of Peptide Vaccine Design). Let V = {v1, v2, · · · , vn} and 2 ≤ k ≤ n. The similarity graph G(V, E) is shown in Figure 1 where the two nodes connected by each edge are considered similar. Each node in the graph has a degree of at most 1, i.e., it can form only one edge, except for v1 which has a degree of 2. The objective function f has properties as follows:

• For any feasible solution X ⊆ V with |X| ≤ k - 1, and ∀vx ∈ V \ X such that ∀vy ∈ X, (vx, vy) /∈ E, it holds f (X ∪ {vx}) = f (X) + f ({vx}); (4)
• For any vx, vy ∈ V \ {v1, v2, v3}, v ∈ {v2, v3}, it holds f ({vx}) = f ({vy}) < f ({v}) < f ({v1}); (5)

Theorem 2. For the peptide vaccine design example given in Definition 3, PVD-GSEMO-WR can find an optimal solution within an expected number of iterations O(kn2 ), while the greedy algorithm Optivax-P cannot.

Proof. We first analyze the property of the optimal solutions, and then show that Optivax-P cannot find an optimal solution. Finally, we derive the expected number of iterations required for PVD-GSEMO-WR to find an optimal solution.

Let Xi ⊆ V denote a feasible solution with |Xi | = i and the set {v2, v3} ⊈ Xi . For any X2 and vx ∈ V such that (v1, vx) /∈ E, it holds that

where the equalities hold by applying Eq. (4), the first inequality holds by Eq. (5), that is, ∀v ∈ V, f({v}) ≤ f({v1}), and the last inequality is by Eq. (6). This implies that the solution {v2, v3} is the unique optimal solution for k = 2. When 2 < k < n, we can similarly derive that for any Xk =X2 ∪Xk-2 , it holds that

where the equality holds by Eq. (4), and the inequality holds by Eq. (7). Here Xk-2 is assumed to not contain the peptides v1, v2 and v3. If v1 ∈ Xk-2 , we can exchange v1 with a peptide v fromX2 to ensure v1 /∈ Xk-2 . The same exchange process can be applied if v2 ∈ Xk-2 or v3 ∈ Xk-2 . Note that by the definition of Xi , we know that Xk can contain at most one of the three peptides v1, v2 and v3; otherwise, Xk will violate the pairwise constraints in Figure 1 or the condition {v2, v3} ⊈ Xi . Eq. (5) states that for any vx, vy ∈ V \ {v1, v2, v3}, f({vx}) = f({vy}); thus, f(Xk-2 ) has a constant value, denoted as C ≥ 0. Combining Eq. (8), we obtain that f(Xk ) < f({v2, v3})+C. Thus, we can find that a feasible solution O is optimal iff O = {v2, v3} ∪ Xk-2 .

The greedy algorithm Optivax-P first selects v1, which has the largest marginal gain. Next, Optivax-P deletes all peptides from the set V that are connected to v1, including v2 and v3. Thus, the optimal solution cannot be found.

For the PVD-GSEMO-WR algorithm, the problem is implemented as maximizing f(X) and minimizing |X| simultaneously. We then prove that PVD-GSEMO-WR can find an optimal solution O = {v2, v3} ∪ Xk-2 benefiting from the bit-wise mutation operator and repair strategy, respectively. Let Oi (2 ≤ i ≤ k) denote the best feasible solution with |Oi | = i. We can verify that Oi must contain v2 and v3, and inserting a specific element into Oi can generate Oi+1 . We reiterate that the definition of Xi ⊆ V represents a feasible solution with |Xi | = i (0 ≤ i ≤ k) and the set {v2, v3} ⊈ Xi . In the following proof, we first show that how PVD-GSEMO-WR can find Oj from Xi (i < j ≤ k), and then follow the path Oj → Oj+1 → · · · → Ok = O to produce an optimal solution.

[Bit-wise mutation] When Xi does not contain the peptide v1, the idea is that flipping at most two 0-bits (corresponding to v2 and v3) of the solution Xi and keeping other bits unchanged can find Oi+1 or Oi+2 . Note that Xi may contain v2 or v3, but will not contain them simultaneously due to the condition {v2, v3} ⊈ Xi . As the solutions in the population P are incomparable and f2(X) = -|X|, P contains at most one solution for each subset size 0, 1, · · · , k. Furthermore, the empty solution X0 is guaranteed to be present in the initial population. Thus, the probability of generating Oi+1 or Oi+2 from Xi is at least (1/|P|) · (1/n)2 · (1 - 1/n)n-2 ≥ 1/(en2 |P|), where 1/|P| is the probability of selecting Xi by uniform selection in line 4 of Algorithm 2, and (1/n)2 · (1 - 1/n)n-2 is the probability of flipping two specific 0-bits (corresponding to v2 and v3) of the solution Xi while keeping other bits unchanged in line 5. Because |P| ≤ k + 1, the expected number of iterations for finding Oi+1 or Oi+2 is at most en2 |P | ≤ e(k + 1)n2 = O(kn2 ).

[Repair strategy] We now show how PVD-GSEMO-WR escapes a local optimal solution Xi (1 ≤ i < k) that includes v1. This is achieved by employing the repair strategy. The process begins by simultaneously flipping the two 0-bits in Xi corresponding to v2 and v3. This action generates an infeasible solution v1 ∪ Oi+1 , which now contains the edges (v1, v2) and (v1, v3) in Figure 1. Subsequently, the repair strategy, as specified in lines 3-4 of Algorithm 4, randomly chooses to remove either v1 or v2. Assuming it opts to exclude v1, we then obtain the solution Oi+1 . Thus, the probability of generating Oi+1 from Xi is at least (1/|P |) · (1/n)2 · (1-1/n)n-2 · (1/2) ≥ 1/(2en2 |P |), where 1/2 is the probability of excluding v1 from {v1, v2} in line 5 of Algorithm 4. This implies that the expected number of iterations is at most 2en2 |P | = O(kn2 ).

Taking the maximum of the expected number of iterations, Oj can be generated from Xi (i < j) within O(kn2 ) expected number of iterations. Note that once Oj (2 ≤ j ≤ k) is generated, it will always be kept in P, since it cannot be dominated by any other solution. The probability of the event Oj → Oj+1 is at least (1/|P|) · (1/n) · (1 - 1/n)n-1 ≥ 1/(en|P|), since it is sufficient to select Oj in line 4 of Algorithm 2, and then flip only one specific 0-bit corresponding to a peptide that has no edges with any peptide belonging to the subset Oj in Figure 1. Because the length of the path Oj → Oj+1 → · · · → Ok = O is at most k - 2, the total expected number of iterations for finding an optimal solution O is at most O(kn2 ) + (k - 2) · en|P | = O(kn2 ), implying that the theorem holds.

From the proofs, we can find that the repair strategy improves the ability to jump out of the local optimum. Without the repair strategy, the probability of generating Oj from the local optimal solution Xi containing v1 is at least (1/|P|) · (1/n)3 (1 - 1/n)n-3 ≥ 1/(en3 (k + 1)), where (1/n)3 (1 - 1/n)n-3 is the probability of flipping two 0-bits (corresponding to v2 and v3) and one 1-bit (corresponding to v1) of Xi while keeping other bits unchanged. Thus, the expected number of iterations to escape from local optima is O(kn3 ), which is only O(kn2 ) if using the repair strategy.

5 Empirical Study

In this section, we examine the performance of PVD-EMO on a peptide vaccine design for COVID-19, by comparing its two variants, PVD-GSEMO-WR and PVD-NSGA-II-WR, with Optivax-P, the state-of-the-art greedy algorithm that has outperformed thirty other algorithms in the peptide vaccine design as reported in [Dai and Gifford, 2023]. Note that PVDGSEMO-WR and PVD-NSGA-II-WR correspond to PVDEMO using GSEMO and NSGA-II, respectively.

The experiments are mainly to answer two questions: Whether any variant of PVD-EMO is better than the previous algorithm Optivax-P? Can the warm-start and repair strategies improve efficiency and performance, respectively?

The initial populations of the two variants of PVD-EMO are generated using the warm-start strategy. Inspired by the recent theoretical work [Zheng et al., 2022], the population size of PVD-NSGA-II-WR is set to 2 times the size of the Pareto front, i.e., 2(k + 1). The warm-start strategy creates one feasible solution for each size 0, 1, · · · , k for PVDGSEMO-WR, and two feasible solutions for each size for PVD-NSGA-II-WR, where the output of Optivax-P is used as one of the solutions for size k. PVD-NSGA-II-WR applies one-point crossover and bit-wise mutation with probabilities of 0.9 and 1, respectively. As PVD-EMO is an anytime algorithm, whose performance will be gradually improved by increasing the number of iterations, we set the number of objective evaluations to 20kn, to make a trade-off between the performance and runtime, compared to kn evaluations used by Optivax-P. Also it is randomized, and thus we run it ten times independently and report the average values.

We use the same dataset of producing a peptide vaccine for COVID-19 as in [Dai and Gifford, 2023], which consists of a set of candidate peptides (|V | = 1043), a set of genotypes (|M| = 1018459) for Major Histocompatibility Complex class I (MHC-I), their frequencies w(m) derived from diverse populations to be representative of the global population, and the binding probability pv,m for each peptide-MHC pair, generated by the SOTA neural network-based model NetMHCpan [Reynisson et al., 2020]. The set of pairwise constraints require that no pair of peptides can be within 6 edits (insertions, deletions, or substitutions) of each other.

We set k ∈ {30, 40, . . . , 70} and N = ⌊0.25k⌋. The results are shown in Figure 2. It can be observed that the performance of PVD-EMO will be at least as good as the competitive greedy algorithm Optivax-P. For a relatively simple problem (k = 30), Optivax-P may find the optimal solution, in which case PVD-EMO performs equally well. However, as the problem complexity increases (k ∈ {40, 50, 60, 70}), both PVD-GSEMO-WR and PVD-NSGAII-WR surpass Optivax-P, showing the superiority of the PVD-EMO framework. This may be because PVD-EMO naturally maintains a population of diverse solutions due to the bi-objective transformation, the employed bit-wise mutation operator has a good global search ability, and the repair strategy further improves the search ability. These characteristics can lead to a better ability of escaping from local optima. Among the two variants of PVD-EMO, PVD-GSEMOWR performs better than PVD-NSGA-II-WR, which may be because the population of NSGA-II may contain redundant dominated solutions, leading to the bad performance.

To more clearly examine the effectiveness of bi-objective reformulation, we conduct experiments using the singleobjective EA (µ + 1)-EA. For fair comparison, (µ + 1)- Runtime in kn EA, with a population size of µ = k + 1, follows the settings of PVD-GSEMO-WR and is labeled as (µ + 1)-EAWR. The key difference is that (µ + 1)-EA-WR updates its population by retaining the best µ solutions according to their objective f values, whereas PVD-GSEMO-WR preserves non-dominated solutions due to the bi-objective reformulation. Additionally, to assess the impact of warm-start and repair strategies, we test two variants: PVD-GSEMO-R, which solely employs the repair strategy, and PVD-GSEMO with no strategies. We plot the curve of the objective value over runtime with k = 40 in Figure 3. Optivax-P is a fixedtime (nearly kn) algorithm, while the others are anytime algorithms that can achieve better performance with increased runtime. The results show that even without warm-start, both PVD-GSEMO-R and PVD-GSEMO outperform Optivax-P within 4kn and 7kn, respectively. PVD-GSEMO-WR and PVD-GSEMO-R converge to the best objective value, while PVD-GSEMO-WR is faster due to the warm-start advantage. In comparison, (µ + 1)-EA-WR, despite employing both strategies, attains the second-best objective value, thus validating the value of bi-objective reformulation to maintain a diverse population. Meanwhile, PVD-GSEMO, lacking the repair strategy, only achieves the third-best objective value, underscoring the vital role of the repair strategy in avoiding local optima, as demonstrated in Theorem 2. These findings validate the effectiveness of bi-objective reformulation along with the warm-start and repair strategies.

Figure 3: The average objective value ± the standard deviation vs. runtime (i.e., number of objective evaluations) with k = 40.

6 Conclusion

This paper proposes the PVD-EMO framework, employing any MOEA to solve the bi-objective reformulated PVD problem. We incorporate warm-start and repair strategies to improve efficiency and performance. We prove that the warm-start strategy ensures that PVD-EMO maintains the same worst-case approximation guarantee as the state-of-theart greedy algorithm Optivax-P. Furthermore, we prove that unlike Optivax-P, PVD-EMO can successfully avoid getting stuck in local optimal solutions. Empirical results on a peptide vaccine design for COVID-19 show that PVD-EMO can achieve better performance than Optivax-P. An interesting future work is to design better MOEAs (e.g., using fast mutation operator [Doerr et al., 2017]) for PVD-EMO.

Acknowledgments

The authors want to thank the anonymous reviewers for their helpful comments and suggestions. This work was supported by the National Science and Technology Major Project (2022ZD0116600) and National Science Foundation of China (62276124). Chao Qian is the corresponding author.