A Minimal Framework for Optimizing Vaccination Protocols Targeting Highly Mutable Pathogens
Saeed Mahdisoltani Pranav Murugan Arup K Chakraborty Mehran Kardar









saeedmah@mit.edu
Department of Chemical Engineering, Massachusetts Institute of Technology, Cambridge, MA 02139
Institute for Medical Engineering and Science, Massachusetts Institute of Technology, Cambridge, MA 02139
Department of Physics, Massachusetts Institute of Technology, Cambridge, MA 02139


Department of Physics, Massachusetts Institute of Technology, Cambridge, MA 02139
Department of Electrical Engineering and Computer Science, Massachusetts Institute of Technology, Cambridge, MA 02139


arupc@mit.edu
Department of Chemical Engineering, Massachusetts Institute of Technology, Cambridge, MA 02139
Institute for Medical Engineering and Science, Massachusetts Institute of Technology, Cambridge, MA 02139
Department of Physics, Massachusetts Institute of Technology, Cambridge, MA 02139
Ragon Institute of Massachusetts General Hospital, Massachusetts Institute of Technology and Harvard University, Cambridge, MA 02139
Department of Chemistry, Massachusetts Institute of Technology, Cambridge, MA 02139


kardar@mit.edu
Department of Physics, Massachusetts Institute of Technology, Cambridge, MA 02139







    A persistent public health challenge is finding immunization schemes that are effective in combating highly mutable pathogens such as HIV and influenza viruses. To address this, we analyze a simplified model of affinity maturation, the Darwinian evolutionary process B cells undergo during immunization. The vaccination protocol dictates selection forces that steer affinity maturation to generate antibodies. We focus on determining the optimal selection forces exerted by a generic time-dependent vaccination protocol to maximize production of broadly neutralizing antibodies (bnAbs) that can protect against a broad spectrum of pathogen strains. The model lends itself to a path integral representation and operator approximations within a mean-field limit, providing guiding principles for optimizing time-dependent vaccine-induced selection forces to enhance bnAb generation. We compare our analytical mean-field results with the outcomes of stochastic simulations and discuss their similarities and differences.






A Minimal Framework for Optimizing
Vaccination Protocols 
 Targeting Highly Mutable Pathogens
    Mehran Kardar
    April 29, 2024
=============================================================================================







§ INTRODUCTION
	
    

The adaptive immune system has a remarkable ability to detect and combat a virtually unlimited number of previously unseen pathogens <cit.>. Immune responses are mediated by T lymphocytes (T cells) and B lymphocytes (B cells), each equipped with surface receptors known as T cell receptors (TCRs) and B cell receptors (BCRs), respectively. 
In a healthy human adult, there are approximately 10^11 T cells and B cells distributed throughout the body <cit.>. 
With the staggering number of TCR and BCR sequences (> 10^14) generated via the V(D)J recombination process <cit.>, most T cells and B cells express a unique surface receptor distinct from others. This extensive receptor diversity enables the immune system to mount specific responses against each new pathogen encountered, including those that emerge after an individual's birth. 


Upon exposure to a new pathogen or vaccine component (collectively referred to as antigens), B cells and antibodies with high binding affinity for the surface proteins of the pathogen are generated by a process called affinity maturation <cit.>. 
Affinity maturation entails an accelerated Darwinian evolution of naive B cells, and it occurs within specialized structures called germinal centers (GCs), which are transiently formed in secondary lymphoid organs such as lymph nodes <cit.>. 
First, naive B cells undergo activation and become a candidate for entry into GCs if their BCRs can bind sufficiently strongly to specific residues known as epitopes on the surface proteins of the antigen <cit.>. 
Upon entering the GC, the activated B cells undergo rapid replication while concurrently accumulating somatic mutations into their BCRs at an accelerated rate. 
The GC B cells then interact with the infecting antigen, which is displayed on Follicular Dendritic Cells (FDCs) located within the GC environment.  
B cells whose receptors exhibit stronger binding to the displayed antigen have a greater probability of internalizing the antigen and are consequently more likely to undergo positive selection. Positive selection involves multiple stages of productive interactions of the B cells with certain T cell types within the GC. 
Conversely, B cells with low-affinity BCRs are typically eliminated via apoptosis as they are less likely to successfully compete with high-affinity B cells for the limited selection signals <cit.>. 



The majority of positively selected GC B cells undergo multiple rounds of replication, somatic mutation, and selection processes. As a result of the competition among GC B cells during the affinity maturation process, an initially naive B cell population with low-affinity BCRs gradually evolves to exhibit strong binding to the antigen. 
During each round, a small number of the positively selected B cells exit the GCs and differentiate into either memory B cells or plasma cells <cit.>. 
Plasma cells are responsible for secreting antibodies, which are soluble forms of the BCR that can bind the antigen's surface proteins with high affinity and neutralize the pathogen's ability to infect host cells. Thus, the affinity maturation process generates an effective immune response by producing antibodies that specifically neutralize a particular pathogen <cit.>.


Memory B cells, on the other hand, play a crucial role in the immune system by retaining a memory of the antigen that activated their parent B cell during the initial exposure. This memory enables a rapid immune response upon re-encountering the same antigen in the future <cit.>. 
Upon re-exposure to a pathogen, existing memory B cells are selected in an affinity-dependent manner and rapidly expanded outside GCs, a process akin to the mechanisms occurring within GCs but with little to no mutation <cit.>. A significant fraction of these expanded memory B cells differentiate into plasma cells,  secreting a surge of antibodies that offer immediate protection. Meanwhile, new GCs are also formed, producing even higher affinity memory B cells and antibodies over more extended time scales, thus bolstering the immune response against the recurring pathogen.  


Vaccines leverage immunological memory to confer protection against severe or life-threatening illnesses caused by specific pathogens they target <cit.>.
By eliciting memory cells and antibodies, vaccines prime the immune system to mount a rapid and specific response to the targeted pathogen upon subsequent exposures <cit.>. 
While vaccines have significantly improved in terms of safety and efficacy over time, developing effective vaccines against highly mutable pathogens remains a continuing challenge. Pathogens such as HIV, influenza, and SARS-CoV-2 mutate rapidly, allowing them to evade the immune responses targeted at specific strains.  

Surface proteins of highly mutable viruses can exhibit significant variation across different strains, rendering vaccine-induced antibodies that neutralize a specific strain ineffective against other strains. This variability complicates vaccine development efforts, as vaccines must contend with the ever-evolving nature of these pathogens.

Certain residues within the surface proteins of highly mutable viruses must remain relatively conserved across mutant strains to preserve the virus's ability to infect host cells. 
Antibodies targeting these conserved residues hold the potential to neutralize a broad spectrum of virus strains <cit.>. 
In some HIV-infected patients, broadly neutralizing antibodies (bnAbs) naturally evolve, demonstrating the immune system's capability to elicit such cross-reactive antibodies <cit.>. 
However, bnAbs are seldom generated in response to natural infection or conventional vaccination, and their emergence typically occurs over extended periods and in limited quantities <cit.>. 
Vaccination strategies aimed at eliciting bnAbs have therefore emerged as a promising approach for diseases like HIV and influenza <cit.>. 
These vaccines have the potential to confer protection against variant strains by inducing a robust immune response targeted at conserved epitopes. 
Despite the advancements, designing universal variant-proof vaccines for mutating pathogens remains a formidable challenge, as the typical immune response to such pathogens is predominantly driven by strain-specific B cells and antibodies.   


For modeling purposes, we can conceptualize the affinity maturation process as the evolutionary population dynamics of B cells `species' within GCs, characterized by inherent stochasticity stemming from randomness in various steps including B cell activation, GC entry, replication, antigen internalization, T cell selection, and more.  
In this framework, vaccine antigens serve as external interventions  guiding this stochastic dynamics process towards desired goals. Designing vaccination protocols that confer broad protection against mutant viruses entails selecting appropriate external factors that steer GC processes to increase the likelihood of rare evolutionary trajectories leading to the development of bnAbs. For instance, external manipulations can involve immunization with variant antigens sharing conserved regions of a virus' spike while differing in variable parts. Experimental and computational investigations have focused on enhancing bnAb production through carefully designed vaccination protocols <cit.>, but achieving success remains elusive. Albeit constrained, there are various routes to guiding affinity maturation, including use of variant antigens, number of vaccine doses, timing between shots, antigen composition in each shot, and dosage. The large space of variables for choosing vaccination protocols makes an unguided search unlikely to succeed. 

Computational approaches have proven invaluable in systematically exploring the large space of vaccination protocols, providing insights into the outcomes of affinity maturation in response to diverse vaccine designs <cit.>. 
Several computational studies have investigated the efficacy of sequential and cocktail immunization procedures for eliciting bnAbs, highlighting an optimal difference between  the variable parts of antigen variants included in vaccine shots that maximize the likelihood of bnAb production <cit.>. 
Signatures of such optimal antigen differences have been observed in HIV patients who naturally develop bnAbs <cit.> and in mice vaccinated by sequentially administered influenza antigens <cit.>. 
Mechanistic explanations for why an optimal difference between variant antigens promotes bnAb evolution have also been proposed <cit.>. 
However, existing studies typically focus on immunization with a fixed  cocktail of antigens <cit.> or  sequential procedures where the variable parts of the antigens are substantially different, maximizing the selection pressure in favor of bnAbs <cit.>. 





In this paper, we use analytical and computational methods to model the complex task of designing vaccines capable of eliciting immune responses that confer protection against diverse mutant strains within a virus family. 

To unveil guiding principles for optimizing bnAb production, we investigate generic time-dependent vaccination schemes that involve administering mixtures of relatively similar antigens over time. 
Our model significantly expands upon previous works in Refs. <cit.> and <cit.> by allowing variations in both the location and spread of the fitness landscape induced by the vaccination protocol. To describe BCR-antigen interactions, we adopt a shape-space representation <cit.>, which enables us to capture the essence of the evolutionary dynamics of B cell populations in response to vaccination in a minimal framework. 



The rest of the paper is organized as follows: In Section <ref>, we describe a simplified model that captures the essential steps of affinity maturation of the B cell population in GCs, namely replication, apoptosis, and mutations. After developing a mean-field approximation to the stochastic dynamics, Sections <ref> and <ref> provide optimization results based on a path integral representation, and operator approximation methods. In particular, for a Gaussian initial germline B cell population that is away from the bnAb state, we derive an exact identity (Eq. (<ref>)) that directly provides the mean-field prediction for the optimal location for the focus of the fitness landscape imposed by the vaccine (Eq. (<ref>)). Using this identity, along with operator methods for approximating the solution to the population dynamics, we obtain a result (Eq. (<ref>)) that enables us to investigate the optimal spread of the vaccine fitness. We next 

compare the mean-field predictions with stochastic simulations in Section <ref>. Finally, we present concluding remarks and outline future research directions in Section <ref>. 
There are four appendices that contain further information regarding the derivation of the continuum mean-field model (Appendix <ref>), the details of the path integral calculation (Appendix <ref>), the operator approximation for the discrete-bin dynamics (Appendix <ref>), and different aspect of the B cell dynamics under different vaccine protocols (Appendix <ref>).  



§ THEORETICAL MODEL
 
    






 §.§ Minimal model for BCR evolution
 


To construct a mechanistic model of affinity maturation, we frame it as the population dynamics governing the evolution of B cells–or equivalently, their BCRs–where replication rates are determined by their binding affinity to the stimulating antigen. 
The interaction between BCRs and antigens is typically captured through phenomenological approaches such as the string model <cit.>. 
Here we adopt a more coarse-grained approached based on the concept of shape-space representation <cit.>.  

In particular, BCRs and antigens are represented by points in a  d-dimensional Euclidean space (Fig. <ref>). Each dimension corresponds to parameters pertinent to calculating BCR-antigen interactions, including factors such as amino acids sequences, spatial conformations, hydrophobicity, and the like. 
As these coordinates in the shape space can assume distinct values, we initially discretize the shape space into a number of similarity bins indexed by x=(x_1,x_2,…,x_d). 

BCRs and antigens positioned closely in shape space are presumed to have complementary characteristics that allow them to bind strongly, while those situated farther apart exhibit weaker binding <cit.>. 

Following the approach outlined in Ref. <cit.>, we characterize the binding free energy between a BCR and an antigen located in similarity bins R_BCR and R_Ag by their distance as  

    E_bind∝ k_BT    ||R_BCR -R_Ag||^2.
 
Consequently, the associated binding affinity is expressed as (E_thr-E_bind), with E_thr denoting the threshold binding free energy for B cell activation. 
B cells bind to the antigen presented on the FDCs with an affinity determined by the equilibrium constant of the BCR-antigen binding, defined as  
K_a = exp[  (E_thr-E_bind)/k_BT]. 
Increased affinities correspond to  stronger BCR-antigen binding with higher equilibrium constants. Affinity maturation can amplify the binding affinity of naive BCRs by tenfold and its associated K_a by a thousandfold <cit.>. 
Furthermore, we assume that BCRs located closer to the origin (x=0) in the shape space have greater affinity toward the conserved residues of the antigen surface proteins. 
These BCRs exhibit increased  tolerance to mutations in the surrounding variable residues of the antigen, facilitating the elicitation of bnAbs. Our objective, therefore, is to devise a vaccination procedure conducive to directing affinity maturation to optimally guide the initially naive BCRs toward the x=0 bin.

Since high-affinity B cells are more likely to undergo positive selection and further replication, the binding affinity to the vaccine antigen(s) effectively imposes a fitness function on GC B cells. 
When a single antigen type is included in a vaccine shot, the induced fitness function exhibits a sharply peaked distribution centered around the complementary BCR bin within the shape space. This distribution features a minimal spread, denoted as σ_min, to encompass nearby bins whose BCRs may also bind the given antigen, albeit with reduced affinity. 
Conversely, in the case of a cocktail shot containing mixtures of different antigen types, the overall fitness function results from the summation of the individual fitness functions induced by each antigen type.[We are assuming a `see all' scenario wherein all antigens are well-mixed and presented homogeneously on FDCs, such that B cells encounter all the vaccine antigens at every round <cit.>.] 
This typically yields a non-convex fitness landscape with multiple peaks and valleys. 
However, if the antigens included in a cocktail shot are sufficiently similar (e.g.,  sharing a common epitope containing conserved residues and surrounding variable portions that are different), the individual fitness peaks would cluster closely in shape space. 
Henceforth, we assume the administered antigens at each time point are sufficiently similar to allow for approximating the overall vaccine-induce fitness landscape with a smooth convex function encompassing these peaks (Fig. <ref>). The resultant coarse-grained fitness function is then roughly centered in the midst of the individual antigens, with a spread generally narrower for more similar antigens. 
We represent this coarse-grained fitness function as V(x), with its peak located at x_V (referred to as the vaccine center or focus) and its spread denoted as σ_V. 
It is worth noting that under the assumption of similar antigens, x_V may generally be distant from the BCR bin associated with the bnAb sequences. 







Now let us describe a minimal model for the affinity maturation process within the representation outlined above. 
The germline distribution that initially seeds the GCs is denoted by n(x,t=0) ≡ n_0(x). As affinity maturation progresses, this distribution evolves over time to n(x,t). We model the evolution of BCRs during affinity maturation through the following stochastic processes: 

    
  * Replication:  the rate of B cell replication is determined by the stimulating antigen(s) included in the vaccine dose. By altering the composition of the presented antigens, we assume it is possible to modulate the vaccine-induced fitness, V(x,t), across different bins in  shape space, with corresponding time-dependent vaccine center and spread,  x_V(t) and σ_V(t).  
    Note that we introduce a time dependence in V to accommodate vaccination protocols where different antigens or cocktails of antigens are introduced at different times.  
    
  * Mutation: BCR mutations occur with rate γ_x,y, corresponding to jumps between shape-space bins x and y. 
    While nucleotide mutations can theoretically move a BCR between any two similarity bins, the coarse-grained nature of shape-space coordinates–reflecting the effect of many amino acids on the binding affinity–suggests that large jumps are exceedingly rare.
    For the simplified mean-field model in the next section, we will assume mutational jumps are local and only occur between nearby bins. This assumption allows us to model them by the diffusion of BCRs in the appropriate continuum limit.   
    
  * Apoptosis:  B cells die with rate λ_x, representing their baseline apoptotic tendency within GCs when not receiving positive selection signals <cit.>. 
  


The rules outlined above describe a generalized linear birth-death-mutation process and are encapsulated in the following Fock space master equation <cit.>

    ∂_t P({ n } , t) = 
             ∑_x V_x(t)  [ 𝔼_x^--1 ] (n_x P({ n } ,t) )
            
      +   ∑_x, yγ_x,y[ 𝔼^+_x𝔼_y^- -1 ] 
        ( n_x  P({ n } , t) ) 
         
      +   ∑_xλ_x [𝔼^+_x - 1] (n_x  P({ n } ,t) )   .


Here, the BCR population is denoted by
{ n } = { n_1, n_2  …} 
where n_i is a non-negative integer representing the occupation number of the ith bin, and P({ n }, t) is the probability of having that specific population at time t. 
On the right-hand side of Eq. (<ref>), we introduce the operators 𝔼_x^±, defined as  

    𝔼_x^±(n_y P({ n } , t) ) ≡ (n_y±δ_x,y)   P({ n}±𝕀_x)   ,
 
where x and y are bin indexes, δ_x,y is the Kronecker delta,  and { n }±𝕀_x denotes a population that differs from { n } by one more or less B cell residing in bin x.  

Each line on the right-hand side of Eq. (<ref>) represents the total rate of change of population probabilities (i.e., gain minus loss) caused by replication, mutation, and apoptosis events, respectively.  
It is worth noting that by describing the occupancy numbers of the shape-space bins, the solution to Eq. (<ref>) would provide a probabilistic description of the entire BCR population in the high-dimensional shape space.   





    


 §.§ Mean-field dynamics and the continuum limit
    


The abstract model described above offers a simplified perspective on affinity maturation. However, deriving insights from the resulting master equation (Eq. (<ref>)) is challenging due to the exponentially large Fock space and the interconnectedness of the time-evolution for all configurations {n}. 

To address this, and obtain an analytically tractable formulation, we consider the dynamics of population averages, denoted by
n_x(t) = ∑_{ n } n_x P({ n },t), and, subsequently, taking the continuum limit of the shape space. 
This leads to the following continuum mean-field equation (see Appendix <ref> for derivation) 

    ∂_t n(x,t) = 
        D  ∇^2 n(x,t)
        + [V(x,t)- λ]   n(x,t)    ,

where, for notational simplicity, we retain the symbol n to represent the average population density in the continuum limit of x.  
We also assume the death rate λ and the diffusion coefficient D (obtained from the underlying mutation rates via Eq. (<ref>)) are constants independent of x.  
Note that, as described earlier,   we consider vaccines with either single antigen types or cocktails wherein the administered antigens are not substantially different from each other. Consequently, at any given time point, the most cross-reactive BCRs to the administered vaccine antigens–namely, BCRs located at x_V(t) in the center of the fitness–may not correspond to a bnAb sequence. This is somewhat distinct from the setting described in Ref. <cit.>, where antigens are different enough that bnAb sequences are always the most fit. 

Solving the mean-field dynamics (Eq. (<ref>)) requires specifying the initial germline distribution, which typically has little overlap with the target bin in  shape space. 
We thus assume n_0(x) is centered at a distant location μ_0 and has a width σ_0 that is an indication of the diversity of the germline BCRs (see Fig. <ref>). 
In the spirit of the central limit theorem, we model n_0(x) by a normal distribution

    n_0(x) = N_0  e^-(x-μ_0)^2/2σ_0^2(2πσ_0^2)^d/2≡ N_0 
        𝒢(x-μ_0, σ_0^2)  ,

where  𝒢 denotes the Gaussian function. 
Here, N_0 represents the total number of activated B cells at t=0.   

Our objective is to determine an optimal vaccination strategy, characterized by  
{x_V(t),σ_V(t)},  
that maximizes the final bnAb population denoted by n(x=0,t_f) ≡ n_f(0), while adhering to a constrained total vaccine dose, i.e.

    ∫_0^t_f dt ∫ d^dx  V(x,t) = constant .


Equation (<ref>) resembles the imaginary-time Schrodinger equation with a time-dependent potential <cit.>. Typically, such equations lack analytical solutions, rendering a direct approach to the maximization problem infeasible.  





Before delving into the analsis of the model, we summarize the assumptions made to simplify the model and acknowledge their associated limitations:  


    
  * We interpret Eq. (<ref>) as imposing a constraint on the total amount of administered antigens in vaccine doses. 
    This interpretation implies a linear relation between the fitness function V (appearing on the left-hand side of Eq. (<ref>)) and the amount of the administered antigen. However, V(x,t) can generally be a complicated function of the antigen concentration, incorporating nonlinear terms that arise from various factors such as competition between B cells for GC entry, positive selection within GCs, and epitope masking (see, e.g., Refs. <cit.>). 
    Nonetheless, the linear relation between fitness and vaccine dose remains a plausible approximation if selection based on the amount of internalized antigen is not overly stringent. This approximation is consistent with the linear structure of the birth-death-mutation model of Eq. (<ref>). 
    
  * We assume that mutational jumps occur independently of replication events, exert incremental effects on BCR affinity, and they are equally probable in any direction (i.e., they are isotropic).    
    These assumptions underlie the symmetric form of the diffusion term in the mean-field dynamics, as included in Eq. (<ref>), along with the vanishing drift term (see Appendix <ref>). 
    In a more realistic model, mutations would primarily occur during replication events and tend to be deleterious than beneficial, resulting in a coarse-grained dynamics with a nonhomogenous diffusion term and a non-vanishing drift term. In addition, deleterious mutations can also significantly reduce the affinity <cit.>, potentially leading to long-ranged steps in the coarse-grained mean-field dynamics.  





    





§ ANALYTICAL RESULTS
  

    


As mentioned earlier, solving the mean-field equation (<ref>) explicitly for a general vaccine fitness function, and thus direct optimization of bnAbs, is not feasible. In this section, we propose a method to simplify this maximization by applying a coordinate change to the path integral representation of the solution to the mean-field dynamics (Eqs. (<ref>) and (<ref>) below). 







 §.§ Path integral representation and the optimal location of the fitness peak
    




We begin  constructing the path integral formulation by rearranging Eq. (<ref>) for an infinitesimal timestep δ t: 

    n(x,t+δ t)    = { 1 + δ t   D ∇^2 + δ t[V(x,t)-λ] } n(x,t) 
       =  e^δ t D∇^2 e^δ t[V(x,t)-λ] n(x,t) + O(δ t^2)   .
 
Equation (<ref>) provides a formal method for evolving the population profile in infinitesimal timesteps. 
By repeatedly applying Eq. (<ref>) to evolve the population over time, the final population n(x,t_f) ≡ n_f(x) can be expressed by the formal time-ordered product <cit.>


    n_f (x)= lim_M →∞∏_j=0^j=M 
         e^δ t D∇^2  e^δ t [V(x,jδ t)-λ]  n_0(x) .


Here, we have defined δ t = t_f/M with M as the number of time increments, and V_j(x) denotes the vaccine-induced fitness function at time jδ t. 
The products in Eq. (<ref>) act in the operator sense:[This product formula is also known as the first-order Suzuki–Trotter decomposition <cit.> and it forms the basis for the path integral formulation of quantum mechanics <cit.>, quantum simulators <cit.>, and the so-called operator splitting methods <cit.>.] 
starting from the germline population n_0, at the jth infinitesimal timestep δ t, the B cell population is locally expanded/contracted according to the vaccine operator e^δ t [V(x,jδ t) - λ], 
followed by diffusion broadening under the operator 
e^δ t D ∇^2.  

The effect of the diffusion operator on a test function η(x) is given by the convolution 

    e^δ t D ∇^2η(x) = ∫ d^dx' 𝒢(x-x', 2Dδ t)  η(x').


Due to the non-local nature of this term, the final population at any location in  shape space results from all potential ways that different parts of the initial population spread and expand until they reach that location at the final time. 
Upon inserting the diffusion operators in Eq. (<ref>) and integrating over all intermediate points, the final population can be more explicitly expressed as 
   
    ∫∏_j=0^M d^dx_j  𝒢(x_j+1-x_j, 2Dδ t)   
        e^δ t (V(x_j, jδ t)-λ) 
        n_0(x_0) ,

with x_M+1 = x.  
This expression is then rearranged into the standard path integral form <cit.> 


    n_f(x) = e^-λ t_f∫ d^dx_0   n_0(x_0) 
        ∫_r(0)=x_0^r(t_f)=x𝒟r(u)   
        e^-S[r,ṙ] ,


where 
𝒟r∝ dx_1 … dx_M stands for the measure of all paths connecting the endpoints, each weighted according to its `action' 

    S[r,ṙ] = ∫_0^t_f du  ( ṙ^  2(u)/4D - V(r(u),u) ) .
 
Note that Eq. (<ref>) indicates that evolutionary paths passing through regions with high vaccine fitness tend to have larger weights in the path integral expression.  

Equation (<ref>) presents a formal expression for the solution to the mean-field dynamics, applicable to any arbitrary initial distribution and fitness landscape. For the optimization of bnAbs, we substitute the germline distribution (<ref>) into Eq. (<ref>)  and then apply the linear coordinate transformation    
x→x + x^*(t), where


    x^*(t) ≡μ_0/1+σ_0^2/2Dt_f(1-t/t_f)  .


After some algebra as described in Appendix <ref>, and upon comparing the resulting expression with the original one, we derive the following exact relationship


    n_f(0) |_μ_0 = ℓ ; Vn_f(0) |_μ_0 = 0;V 
        = 
        exp - ℓ^2/2(σ_0^2 + 2Dt_f)  .

In this equation, the numerator on the left-hand side denotes the final population at x=0 (i.e., the number of final bnAbs) given that the normally-distributed initial population (<ref>) is centered at μ_0 =ℓ and evolves under the vaccine fitness V(x,t).  
The denominator, on the other hand, represents the final bnAb count for a germline population that is centered at μ_0 = 0 and is evolved under the transformed fitness V defined as 

    V(x,t) ≡ V ( x+x^*(t),t )   .

Therefore, the exact relationship (<ref>) simply relates the number of bnAbs generated by fitness V, with the germline population initially centered away from the target in  shape space, to the number of bnAbs developed under the modified fitness V, assuming the germline is already centered at the target (see top panel of Fig. <ref>).


The crucial observation here is that the right-hand side of Eq. (<ref>) is independent of the choice of the vaccine parameters. Thus, the maximum of the numerator on the left-hand side as V is varied (by changing {x_V(t),σ_V(t) }), must coincide with the maximum of the denominator as V is varied. 
As a result, Eq. (<ref>) effectively converts the maximization of 
n_f(0)|_μ_0=ℓ ; V into the maximization of n_f(0) |_μ_0 = 0; V. 
This mapping simplifies the optimization problem:
since for μ_0 = 0, the starting germline population is already  peaked at, and symmetric around, the target bin  x=0, it makes sense that for V to be optimal, it should also be peaked at the target and symmetric with respect to it at all times (see Fig. <ref>). 
Using this argument, and recalling that the maxima of the numerator and denominator in Eq. (<ref>) coincide, we then conclude that

    x_V(t) |_optimal = x^*(t) 
        =μ_0 2D(t_f-t)/2Dt_f+σ_0^2  .



It is worth noting that for a narrow germline population (σ_0 ≪√(Dt_f)), the optimal vaccine center given by Eq. (<ref>) simply starts close to the naive population at μ_0 and reaches the target state (x=0) at the final time. 
However, for a wider initial population (σ_0 ≳√(D t_f)), the optimal center is shifted closer to the target, reflecting a change in the importance of paths contributing to n_f(0) in Eq. (<ref>).  
Intuitively, one needs to position the vaccine fitness over time such that it gradually moves the germline population from x=μ_0 towards the target state at x=0, while also remaining sufficiently close to the peak of the existing B cell population to avoid  extinction by apoptosis. 
If V(x,t) is positioned too close to the population n(x,t), it may not optimally move the population to reach x=0 within the finite time interval t_f; conversely, if V is centered too far from the peak of n(x,t), there is a risk of population collapse as the  vaccine-induced replication is directed to locations with fewer B cells to begin with. 
The optimal center (Eq.  (<ref>)) aligns with this intuitive understanding and also demonstrates how the germline width σ_0 influences the optimal position of the vaccine fitness. 
In addition, Eq. (<ref>) reveals that the mean-field optimal center moves with a constant velocity from its starting position in the shape space towards the target bin. 
These features are further illustrated in the numerical and simulation results discussed in Section <ref> and Appendix <ref>. 








    

 §.§ Operator approximation and the optimal spread of the fitness function imposed by vaccine antigens
    

As discussed above, Eq. (<ref>) maps the problem of maximizing bnAbs to finding the optimal modified fitness function  
V that maximizes 
n_f(0) |_μ_0 = 0 ; V.  
Upon replacing V by V, Eq. (<ref>)    
describes the expansion of the transformed germline population centered at x=0 under the modified fitness V, coupled with diffusive spreading that arises from BCR mutations. 


In order to investigate the optimal vaccine spread, we employ operator approximations. 
Starting from Eq. (<ref>), we rearrange the operators in the time-ordered product formula  using a second-order truncation of the celebrated Baker-Campbell-Hausdorff  (BCH) formula <cit.>, namely  

    e^ϵ A e^ϵ B = e^ϵ B e^ϵ A e^ϵ^2[A,B](1+O(ϵ^3)),

where A and B are generally non-commuting operators, and C ≡ [A,B] = AB - BA. 
In this formula, we assign 
ϵ→δ t, 
A → D∇^2, and 
B →V_j (x),  
and then utilize this formula to reorder the alternating operator product in Eq. (<ref>) such that all vaccine fitness operators e^δ t  V_j are grouped on the left in the product, all the diffusion operators e^δ t D ∇^2 are collected on the right, and the commutator terms that result from swapping these operators appear in between. The apoptosis operator e^-λ t_f commutes with all the other operators.  


To achieve such an ordering of the operators, we start by moving e^δ t V_0 through the M diffusion operators e^δ t D ∇^2 on its left in Eq. (<ref>), resulting in a net commutator term 
e^(δ t)^2 M C_0, where the commutator operator reads 

    C_0 = [D∇^2 , V_0(x)] = 2 D V_0'(x) ∇ + DV_0”(x)   .

Subsequently, the next fitness operator, e^δ t V_1, needs to pass through (M-1) diffusion operators on its left, thereby creating 
a commutator contribution e^(δ t)^2 (M-1) C_1. 
This process continues for all fitness operators, and in the limit of small increments δ t = t_f/M→ 0, we obtain[Each truncation of the BCH formula introduces an error of order (δ t)^3. It can be shown that the number of such error terms grows at least as ∼ M^3, rendering this truncation an uncontrolled approximation. 
Despite this, the result obtained through this approach matches that of a perturbative solution in powers of the non-dimensionalized diffusion coefficient.   
The second-order truncation remains valid for small diffusion limits and as long as the total time t_f is sufficiently short.]

    n_f  |_μ_0=0;V≈    N_0    e^-λ t_f e^∫_0^t_f dt  V(x,t)
       × e^D∫_0^t_f dt   (t_f - t)  ( 2V'(x,t)∇ + V”(x,t))
       × e^t_f D ∇^2 𝒢(x_0,σ_0^2)   .


Unlike the exact expression in Eq. (<ref>), where the fitness and mutation (diffusion) alternate in the operator product, Eq. (<ref>) suggests that n_f can be approximated by first diffusing the starting profile n_0 under the (full) diffusion operator e^t_f D∇^2. Then, the commutator operators act on this diffused population profile, and, finally, the vaccine fitness operators (and the death operator) expand the profile to its final form. 

The action of e^t_f D∇^2 on the Gaussian initial population  widens its profile and yields 
𝒢(x_0,σ_0^2 + 2Dt_f), which is independent of the fitness profile. 
The symmetry of this updated profile with respect to x_0 = 0 still implies a similar symmetry for the optimal
V(x,t), resulting in V'(x=0,t) = 0. Consequently, Eq. (<ref>) simplifies to  


    n_f(0) |_μ_0=0;V≈      
        N_0 e^-λ t_f 
        e^∫_0^t_fdt V(0,t)
    ×   e^D ∫_0^t_f dt (t_f-t) V”(0,t)  𝒢(0,σ_0^2 + 2Dt_f)   .



For any choice of the vaccine-induced fitness profile, Eq. (<ref>) offers a starting point to investigate the optimal vaccine spread by providing an explicit expression whose optima can be derived by straightforward differentiation.    
Narrowing our focus, we consider the  Gaussian vaccine-induced fitness

    V_𝒢(x,t) ≡ A_V  𝒢(x-x_V(t),σ_V(t)^2 )   .
 
Recall that this choice is consistent with our initial assumption that, in the shape-space representation, the fitness can be parameterized by the location of its center and its spread. Furthermore, Eq. (<ref>) also satisfies the total-dose constraint given by Eq. (<ref>). 
However, it is important to note that Eq. (<ref>) does not represent the most general Gaussian fitness profile allowed by the constraint, as it only considers a fixed vaccine strength A_V.  
The more general case would involve a vaccine dose that varies over time, potentially leading to larger and more diverse antibody responses <cit.>.
We defer the investigation of such scenarios to future studies. 

By selecting the optimal vaccine center according to Eq. (<ref>), the Gaussian fitness of Eq. (<ref>) transforms according to Eq. (<ref>), resulting in the modified fitness
V_𝒢(x,t) = A_V 𝒢(x,σ_V^2(t)), which is still a Gaussian profile but now centered at the origin. 

Substituting V_𝒢 into Eq. (<ref>), we obtain

    n_f(0) |_μ_0=0;V_𝒢≈N_0 e^-λ t_f/√(2π(σ_0^2+2Dt_f))
           ×exp∫_0^t_f dt  ( A_V/√(2π)σ_V(t)(1-D(t_f-t)/σ_V^2(t)) )  .


To maximize this expression, we set its functional derivative with respect to σ_V to zero. Assuming σ_V(t) ≥σ_min at all times (see the discussion below), we finally obtain


    σ_BCH(t) = max{√( 3 D (t_f - t))  , σ_min}


In deriving Eq. (<ref>), we have assumed that σ_V(t) ≥σ_min, which prevents V_𝒢 from reaching arbitrarily large values by imposing a maximum fitness value 
A_V/(√(2π)σ_min).  
As discussed in Section <ref>, the fitness function exhibits a minimal yet finite spread for vaccine shots containing a single antigen type, whereas cocktail shots typically exhibit a wider spread. 
Moreover, in Appendix <ref>, we demonstrate that the optimal spread retains a finite minimum value when employing the BCH approximation for B cell dynamics on discrete similarity bins (before taking the continuum limit in shape space). Additionally, the minimum spread can also reflect practical constraints on how precisely  the vaccine-induced fitness can be concentrated in different regions of  shape space. We posit that σ_min encompasses all of these factors. 

Equation (<ref>) can be understood by noting that within a temporal window of duration Δ t, diffusion affects a range of order √(D Δ t) of bins on the x axis. To maximize n_f(0) while adhering to a constrained total dose, it would thus be sensible to concentrate the fitness on BCRs in a region around the target bin that has the potential to reach the target location by mutation within the available time. 
Replication events at larger mutational distances are unlikely to reach the target within the finite allotted time.
This  suggests that at any time 0<t<t_f, the optimal V should at most cover BCRs within a range of order √(D(t_f - t))  while remaining peaked at x=0. 
Once the vaccine spread reaches the smallest possible value σ_min, it should remain at that value until the final time. 

In a qualitative sense, Eq. (<ref>) suggests that an optimal vaccination protocol should initially have a broader spread and contain shots with more diverse antigens to cover a larger range of BCRs, which would provide the germline population a chance to grow and diversify.  
As time progresses, vaccine epitopes gradually narrow down according to Eq. (<ref>), concentrating the fitness on expanding those BCRs within an accessible mutational distance from the target bin. 


It is also worth noting that the right-hand side of the BCH approximate result (<ref>) depends only on time integrals of V(0,t) and V”(0,t), making it a strictly local function of the transformed vaccine fitness V at x=0. 
This locality aligns with the second-order truncation of the BCH formula, as higher-order terms in BCH's nested commutator expansion would introduce contributions involving more derivatives of the fitness, leading to non-local expressions in x.   

Since the mean-field Eq. (<ref>) aims to describe the coarse-grained temporal evolution of the B cells toward the bnAb state, and considering that bnAbs are typically produced over long timescales in natural circumstances <cit.>, it may be reasonable to assume the slow diffusion regime, where local approximations such as the one above remain reliable at the mean-field level.  
However, it is important to acknowledge that even in the slow diffusion limit, this approximation might break down if either the vaccine fitness or the B cell population profile is sharply peaked, as higher-order derivatives may become significant in the expansion.

Furthermore, as noted in Section <ref>, our approach differs from Refs. <cit.> in our assumption that antigens administered at each time point share sufficient similarity,  that bnAb sequences are not inherently the fittest. 
Given this distinction and the discussion of the previous paragraph, while Eq. (<ref>) bears resemblance to the findings of Ref. <cit.>– where the optimal fitness function appeared narrower in subsequent vaccinations–we approach the interpretation of Eq. (<ref>) with caution.  

In the next section, we will compare the performance of protocols employing a width σ_BCH with those utilizing a narrow fitness width σ_min through numerical results and simulations. The comparative analysis aims to shed light on the effectiveness of different vaccine spread strategies and will elucidate a distinct interpretation of the optimal spread in our study.




§ COMPUTATIONAL RESULTS


     
 




 



In this section, we delve into the dynamics of the B cell population within a one-dimensional shape space (d=1) and investigate both the numerical solution of the continuum mean-field model (Eq. (<ref>)) and simulations of the discrete stochastic dynamics (Eq. (<ref>)). 
It is important to recall that while the master equation deals with discrete shape-space bins, representing BCR occupancy with integer values, the mean-field model offers a continuum description, neglecting population fluctuations. Consequently, it is not immediately clear whether discrete sample averages would align with the mean-field results. 
To establish a comparison between the two approaches, we employ the  stochastic simulation algorithm (SSA) a la Gillespie <cit.> to generate realizations of stochastic dynamics and then compare the sample averages with those obtained from the mean-field treatment. 
This numerical examination sheds light on how effectively the optimal vaccine fitness obtained from the mean-field model translates to the realm of discrete stochastic dynamics


To set up our numerical experiments and simulations, we must define parameters for:  
(1) the germline population, specifically N_0, μ_0, and σ_0; 
(2) the dynamics, including D, t_f, λ; and
(3) the vaccine's strength, denoted by A_V. 
To make these parameter choices more grounded, we draw from empirical observations and rough estimates. 
First, it is cruical to consider the slow and delayed production of bnAbs in HIV patients <cit.>. This implies that the germline population should be initially distant from the target bin (x=0). For instance, it is known that numerous mutations are typically required for germline B cells to evolve into bnAbs that can target CD4 binding sites <cit.>. This informs a constraint  among our parameters: 

    μ_0 ≳√(Dt_f).



Next, in the absence of vaccine-induced replication, the apoptotic B cell population is expected to go extinct before reaching the bnAb state through mutations. In a mean-field approximation, we crudely define extinction as a total population size falling below 1. This extinction condition can be expressed as   

    ln N_0 ≲λμ_0^2/D.
 


Drawing from the findings of Ref. <cit.>, we posit that an effective vaccine would provide sufficient fitness to prevent the extinction.  
We thus choose the vaccine strength A_V so that the total vaccine-induced fitness given at each time exceeds the total death rate across all shape-space bins, i.e., 

    A_V ≳λ N_binΔ x_bin,

where N_bin represents the number of bins and Δ x_bin is the length of a bin, assumed to be one. (Note that as per Eq. (<ref>), A_V has the dimension of (shape-space) length per time.)


Lastly, we note that a further condition to avoid extinction under the optimal vaccine center is for x^* to remain close to the population's center; otherwise, most of the provided fitness is expended on bins with few to no B cells.  
Given the difficulty of analytically solving for the population center as a function of time, we only impose this condition on the germline at t=0, which corresponds to    
|μ_0 - x^*(t=0)| ≲σ_0. 
Substituting for x^*(t=0) from Eq. (<ref>) then yields the fourth condition: 

    t_f ≳σ_0 (μ_0 - σ_0)/D,

which, for a given diffusion constant D, sets a lower bound for the total vaccination time t_f. Typically σ_0 ≪μ_0, ensuring that this bound is significantly smaller than the diffusive timescale μ_0^2/D required for the germline to reach the target by mutations alone (c.f. Eq. <ref>).  

Guided by these conditions, we choose  N_bin = 43 corresponding to the shape-space bins x ∈{ -21,-20,…,0,…,20,21} for the stochastic simulations. This range is associated with the interval (-21,21) for the continuum mean-field numerics.  We set the death rate to  λ = 0.001 and the mutation rate to μ=0.1 in SSA, corresponding to a mean-field diffusion coefficient of D=0.05 (see Appendix <ref>). 
We assume the germline population is normally distributed according to Eq. (<ref>), with an initial population size of N_0=10. In addition, the imposed vaccine fitness follows the Gaussian shape of Eq. (<ref>), where we fix the vaccine strength at A_V=0.08.  
For the stochastic dynamics, we utilize SSA to take averages over 1000 realizations of  Eq. (<ref>) for a total time of t_f = 400. To facilitate the simulation, we implement SSA along with an operator splitting scheme <cit.> where we alternate between periods of stochastic replication events and stochastic mutation or death events (cf. Eq. (<ref>) for the corresponding mean-field expression). In the following, we present the SSA results for M=200 steps, although consistent results are obtained  for M=400. Note that the mean-field numerics presented below are obtained by direct numerical solution of the continuous-time dynamics described given by Eq. (<ref>). 

Fixing the vaccine spread at σ_min=1 and the germline width at σ_0=1, we first examine the effect of different vaccine centers on the final bnAb count by selecting the vaccine center from the family 

    χ^q,p_ℓ,κ(t) ≡ℓ/1+κ^2/2Dt_f(1-(t/t_f)^q )^p   .


In particular, we compare the final bnAb count for χ^1,1_μ_0,σ_0– corresponding to x^* as defined in Eq. (<ref>)–with that of χ^1,1/2_μ_0,σ_0 (vaccine center accelerating towards the target), χ^1/2,1_μ_0,σ_0 (vaccine center decelerating towards the target), and χ^1,1_μ_0,0 (uniform motion similar to x^* but not accounting for germline's width). 
The left panel in Fig. <ref> shows the final bnAb counts as a function of the germline distance from the target for three different vaccine centers x^*(t) (red), χ^1,1/2_μ_0,σ_0(t) (blue), and χ^1/2,1_μ_0,σ_0(t) (green). 
We observe that in each case, the SSA sample averages (solid points) closely match the continuum mean-field numerics (dashed lines). 
To understand this, note that with σ_0 = 1, the germline population is narrowly distributed, so the entire population would roughly move together in  shape space. When this population is relatively far from x=0, moving between neighboring (discrete) bins would only change the distance to the target bin incrementally, thus making the continuum limit a reliable approximation of the discrete dynamics. Additionally, with the chosen vaccine strength (Eq. (<ref>)), the population size increases and makes the population number fluctuations less relevant, thus rendering the mean-field description even more accurate. 
This is also reflected in the fact that the optimal center x^*(t)–obtained from the continuum mean-field approximation–outperforms the other protocols even in the discrete stochastic dynamics. 


The right panel in Fig. <ref> presents a similar comparison of the final bnAb counts, but for a broader germline with σ_0 = 3. Here, we include a fourth choice for the vaccine center, namely χ^1,1_μ_0,0(t) (in yellow), which differs from x^*(t) in that its starting position and velocity are not adjusted by σ_0. We still observe that x^*(t) outperforms the other protocols, although the final bnAb counts (vertical axis) are generally smaller compared to the σ_0 = 1 case. This is because the population does not move toward the target as synchronously when the initial population is spread more widely. 
It is also seen that when μ_0 becomes comparable with σ_0, SSA averages deviate from the continuum mean-field results, which might indicate the breakdown of the continuum approximation for that range of parameters. Nevertheless, the qualitative trend of the SSA sample averages is still captured by the mean-field solution. 

In Appendix <ref>, we present temporal snapshots of the BCR population evolved under various vaccine protocols and explore additional aspects of the population dynamics, such as population size and the location of the population peak over time. Notably, these observations indicate that, as discussed below Eq. (<ref>), the optimal center gradually shifts the germline population from its initial position in  shape space toward the target while maintaining an optimal distance from it. 

Assuming the optimal choice for the vaccine center, we proceed to compare the effectiveness of two different vaccine spread strategies: maintaining the fitness as narrow as possible or adjusting the fitness spread based on the time-dependent BCH result. 
Figure <ref> depicts the final bnAb counts as a function of σ_min. It shows that maintaining the spread at σ_min (in red) results in higher final bnAb numbers compared to setting the vaccine spread according to the time-dependent BCH result given in Eq. (<ref>) (in cyan). 
This observation suggests that maximizing the focus of fitness while aligning its center with the optimal position (as defined in Eq. (<ref>)) results in the largest bnAb numbers ultimately. 
Drawing from the discussion on the minimum fitness spread in Section <ref> (also see Section <ref>), a plausible interpretation of this finding suggests that sequential immunization with a single antigen type, administered at any time, outperforms immunization with cocktails of antigens for which the fitness spread is larger, in agreement with the findings of Ref. <cit.>. This observation implies that even a cocktail that optimally `frustrates' the B cell population, as found in Ref. <cit.>, would be outperformed by a sequential immunization regimen containing appropriately selected single antigen types at each time point. 

It is also noteworthy to contrast Fig. <ref>  with the conclusions drawn in Ref. <cit.>, where the administered vaccine antigens were considered sufficiently distinct in their variable residues such that the fittest BCR sequence corresponded to bnAbs, while the germline would start far from the bnAb sequence. In Ref. <cit.>, the optimal spread of the vaccine-induced fitness function was found to be one that is sufficiently wide to reach the least cross-reactive existing B cells; a narrower spread would have led to B cell extinction. Consequently, the optimal fitness landscape narrowed in subsequent immunizations, a phenomenon that coincides with the narrowing nature of the BCH result (Eq. <ref>). As discussed in the preceding paragraph, the results presented in this work relax the assumption of distinct antigens and instead concern vaccine shots containing similar antigen types.  








 
 




 









 







§ DISCUSSION
    
    We have explored the optimal vaccination protocol for producing bnAbs through a minimal birth-death-mutation model described by the master equation (Eq. (<ref>)) and its continuum mean-field approximation (Eq. (<ref>)).  
Assuming a B cell germline that is initially distant from the desired bnAb state in the shape space, this optimization problem comprises two aspects:  first, determining the optimal time-dependent center of the vaccine-induced fitness function, and, second, its optimal time-dependent spread. 
The center of the vaccine fitness function–i.e., where the fitness peaks–is determined by the epitopes shared among the antigen(s) presented at each time point; the spread of the fitness function, on the other hand,  is controlled by diversity of the presented epitopes in each shot.   
We assumed the administered antigens are sufficiently similar to allow for approximating the vaccine-induced fitness landscape with a coarse-grained smooth function whose peak can be distal from the BCR sequence corresponding to bnAbs (Fig. <ref>). 





We addressed the first part by employing the path integral formulation of the mean-field dynamics (Eq. (<ref>)) and used a linear coordinate transformation to derive a nontrivial exact relationship given in Eq. (<ref>). This relationship is valid for a general vaccine-induced fitness if the germline population is normally distributed (cf. Eq. (<ref>)). 
Equation (<ref>) maps the maximization of the final bnAb count n_f(0) for a germline centered at μ_0 ≠ 0 away from the target, into the simpler problem of finding the optimal transformed fitness, defined in Eq. (<ref>), that maximizes the bnAb number for a germline that is instead initially centered at the target (μ_0=0).  
Based on this mapping, we then argued (see Fig. <ref>) that the optimal vaccine center follows Eq. (<ref>), with x^*(t) defined in Eq. (<ref>). 

The same mapping also allows us to investigate the second part of the bnAb optimization, i.e., the optimal spread of the vaccine-induced fitness landscape. Using a second-order truncation of the celebrated BCH formula, we obtained an approximate form for the mean-field solution (Eq. (<ref>)) that, given a specific form of fitness, can be used to examine how the vaccine spread affects the final bnAb count. In the special case of a Gaussian vaccine-induced fitness (Eq. (<ref>)), maximizing this expression with respect to the vaccine spread led us to Eq. (<ref>). In Eq. (<ref>), we incorporate σ_min as a proxy to encompass all factors constraining the narrowness of the fitness function. These factors may include the (weaker) binding of neighboring BCRs to the presented antigen, the effect of discrete similarity bins, and  experimental limitations. 

Finally, to test the performance of the optimal fitness as obtained from the mean-field model in the discrete stochastic dynamics, we simulated Eq. (<ref>) using the Gillespie SSA and compared the resulting sample averages with the numerical solution of the associated mean-field dynamics. We observed that the SSA averages and mean-field numerics match closely, and the optimal vaccine center, Eq. (<ref>), outperforms other protocols even at the level of the discrete stochastic dynamics (Fig. <ref>). 
In addition, we found that keeping the vaccine spread at the minimum allowed seems to be optimal (Fig. <ref>), suggesting that antigen cocktails reduce the effectiveness of the immunization procedure.   

Our findings suggest that optimal vaccination protocols could potentially be designed to maximize the evolution of rare bnAbs without risking B cell extinction by ensuring that sequentially administered vaccine antigens are not too different at each time point. This important conclusion of our work warrants further investigation by addressing several inherent limitations of the model discussed in Section <ref>, as well as accounting for continuous entry of naive B cells into germinal centers <cit.>.  As briefly discussed in Section <ref>, during affinity maturation, B cells interact with each other in various ways, such as competing for limited survival signals from T cells and memory effects including the epitope masking <cit.>. Incorporating these effects in a birth-death-mutation model  introduces nonlinear terms, and may well lead to modified perspectives on bnAb evolution, compared to the linear model employed in this study. In  future, we intend to explore optimal vaccination strategies using a more realistic version of the stochastic dynamics that accounts for these crucial nonlinearities and feedback loops. Additionally, we plan to extend our analyses to encompass time-dependent vaccine dosages. Ultimately, we hope that our study will motivate further theoretical and experimental investigations by other researchers.





 





§ ACKNOWLEDGMENTS

We would like to thank Leerang Yang and Federica Ferretti for helpful discussions. This research was partially supported by NIH grant 5R01AI175489-02 and the Ragon Institute of MGH, MIT, and Harvard.
MK also acknowledges support from NSF grant DMR-2218849.

Author contributions –
S.M., A.K.C., and M.K. designed the research. S.M. carried out the analytical
calculations, numerical analysis, and stochastic simulations. P.M. provided the original code for Gillespie simulations. S.M., M.K., and A.K.C. analyzed and discussed the results and wrote the paper. PM offered comments on the manuscript.

Declaration of interests –
A.K.C. is a consultant (titled Academic Partner) for Flagship Pioneering and
also serves on the Strategic Oversight Board of its affiliated company, Apriori
Bio, and is a consultant and SAB member of another affiliated company, Metaphore Bio.







§ DERIVATION OF THE MEAN-FIELD IN THE CONTINUUM LIMIT
 
    In this appendix, we elaborate on the derivation of the mean-field equation (<ref>) from the master equation (<ref>). 
The average BCR occupation number at bin x and time t is given by 
n_x(t) = ∑_{ n } n_x  P({ n } , t), where the summation is performed over all possible population configurations { n }.  Note that due to the averaging process,  n_x(t) is not necessarily integer-valued anymore, while x remains discrete at this stage.   
The dynamics of the mean population is obtained by taking a time derivative of n_x(t). Using Eq. (<ref>), it reads

    ∂_t n_x(t) = 
             ∑_{n} n_x∑_x' V_x'(t)  [ 𝔼_x'^--1 ] (n_x' P({ n } ,t) )
        
      +   ∑_{n} n_x∑_x',y'γ_x' y'[ 𝔼^+_x'𝔼_y'^- -1 ] 
        ( n_x'  P({ n } , t) )  
      +   ∑_{n} n_x'∑_x'λ_x' [𝔼^+_x' - 1] (n_x'  P({ n } ,t) )   ,


which can be simplified in each line by an appropriate shift of the variables in the sums over {n}. 
For instance, it is straightforward to show that  

    ∑_{n} n_x∑_x' V_x'(t)  𝔼_x'^- ( n_x' P({n},t) ) = 
           ∑_{n}(n_x + δ_x,x')∑_x' V_x'(t)   n_x' P({n},t)


where δ_x,x' denotes the Kronecker delta. On implementing these simplifications, Eq. (<ref>) turns into the following mean-field dynamics 

    ∂_t n_x(t) =    [V_x(t) -λ_x]  n_x(t)
       +∑_y [γ_yx n_y(t) - γ_xy n_x(t)]   ,

where x and y denote different affinity bins. 

To take the continuum limit of x in Eq. (<ref>), we make the change 

    ∑_y [γ_yxn_y - γ_xyn_x] →
       ∫ d^dΔx [ γ(x-Δx,Δx) n(x-Δx,t) - γ(x,Δx) n(x,t) ]

where γ(x,Δx) denotes the mutation rate for a jump of length Δx that starts from x.  
The assumption of localized jumps corresponds to negligible jump rates for large Δx (see the discussion is Section <ref>), allowing for a Taylor expansion of the integral term in powers of Δx. Truncating the Taylor expansion at second order then gives  

    ∂_t n(x,t)    ≈[ V(x,t) - λ(x) ]   n(x,t) 
        
          + ∇^2 [ D(x)   n(x,t) ] - ∇·[w(x)   n(x,t)],


where we have defined the diffusion coefficient as

    D(x) = 1/2∫ d^dΔx  (Δx)^2  γ(x,Δx)
 
and the drift velocity as w(x) = ∫ d^dΔx Δx γ(x,Δx). 
For symmetric mutation rates, we have w(x,t)=0. Assuming jump rates independent of x then simplifies Eq. (<ref>) into Eq. (<ref>) of the main text. 




§ DETAILS OF THE PATH INTEGRAL CALCULATION
    
     In this appendix, we provide details on the derivation of Eq. (<ref>) from the path integral formula (<ref>). 
We first rewrite Eq. (<ref>) in a new frame, where the origin moves according to
z(t) = z_0 - v t.   
This transformation corresponds adjusting shape-space vectors as
x (t) →x̃ = x-z(t) and paths as 
r(t) →r̃(t) = r(t) - z(t) while keeping the path integral measure unchanged 
(𝒟r̃(t) = 𝒟r(t)). 
In addition, the initial population transforms according to 

    n_0(x_0) = N_0 𝒢(x_0 - μ_0, σ_0^2) 
        
                   →ñ_0(x̃_0) = N_0𝒢(x̃_0 - μ̃_0,σ_0^2)

where we have defined
μ̃_0 ≡μ_0 - z(t=0) = μ_0 - z_0. 





Noting that
ṙ̃̇(t) = ṙ(t) + v, the path integral action (Eq. (<ref>)) transforms as follows

    S[r,ṙ]    = v^2 t_f/4D - v·(r̃(t_f)-r̃(0) )/2D
            +  ∫_0^t_f du (ṙ̃̇^2(u)/4D - V(r̃(u),u) )  
       ≡v^2 t_f/4D - v·(r̃(t_f)-r̃(0) )/2D  +
        S[r̃,ṙ̃̇]

where we have used the fact that 
∫_0^t_f du  ṙ̃̇(u) = r̃(t_f) - r̃_0. 
Here, 
V(r̃(u),u) and S[r̃,ṙ̃̇] denote the transformed vaccine fitness and path integral action which are obtained from V(r(u),u) and S[r,ṙ] by making the change 
r(u) →r̃(u)+z(u). 
Substituting these expressions back into Eq. (<ref>), we obtain for x=0 


    n_f(0) = 
        expv^2 t_f/4D - v·z(t_f)/2D
       × e^-λ t_f∫ d^dx̃_0  ñ_0 (x̃_0)   e^v·x̃_0/2D ∫_r̃(0)=x̃_0^r̃(t_f) = z(t_f) 𝒟r̃(u)   
        e^-S[r̃,ṙ̃̇]  .


To proceed, we utilize the Gaussian shape of ñ_0 and simplify the outer integral by combining the modified initial population ñ_0(x̃_0) = N_0𝒢(x̃_0-μ̃_0, σ_0^2) with the exponential factor, e^v·x̃_0/2D, 
through completing the square, which yields


    ñ_0(x̃_0)   e^v·x̃_0/2D =   N_0𝒢(x̃_0 - μ̃_0 + vσ_0^2/2D , σ_0^2 )   
        e^-v·μ̃_0/2D + v^2 σ_0^2/8D^2  .


We observe that apart from the multiplicative constant,  the transformed path integral expression in Eq. (<ref>) retains the same mathematical structure as the original expression in  Eq. (<ref>). 
 

It is now straightforward to show that for z_0 = v t_f = μ_0/1+σ_0^2/2Dt_f, we have z(t_f)=0 and, moreover,  μ̃_0 - vσ_0^2/2D = 0. Substituting these back into Eq. (<ref>) and using Eq. (<ref>), we finally arrive at


    n_f(0) = exp-μ_0^2/2(σ_0^2+2Dt_f)
         × 
        N_0 e^-λ t_f∫ d^dx_0  𝒢(x_0,σ_0^2) 
        ∫_r(0)=x_0^r(t_f)=0𝒟r(u)   
        e^-S[r,ṙ]  ,


where we have reverted the dummy integration variables x̃_0 and r̃ back to x_0 and r. 
In Eq. (<ref>), the second line on the right-hand side can be recognized as the right-hand side of the original expression in Eq. (<ref>) upon changing V →V and setting μ_0 =0. Rearranging Eq. (<ref>) thus leads to the Eq. (<ref>) in the main text. 

Figure <ref> shows that the left-hand side of Eq. (<ref>), evaluated by numerically solving the one-dimensional mean-field equation for a set of arbitrary vaccine fitness functions, matches the exponential factor on the right-hand side of Eq. (<ref>).  

 








§ BCH APPROXIMATION FOR DYNAMICS ON DISCRETE BINS
 
    Since the SSA are performed on a system of discrete bins, for completeness, this appendix demonstrates how the calculations 
in Section <ref> can be extended to the discrete-bin version of the mean-field dynamics. 

Recall that the general mean-field dynamics for the discretized bins before taking the continuum limit is given by Eq. (<ref>). We set λ_x = λ and γ_xy = γ∑_k=1^d δ_x, y±ê_k where ê_k is the kth unit vector and d denotes the dimension of  shape space. 
Focusing on the one-dimensional case for simplicity, we use the BCH formula (<ref>) with A and B matrices (instead of continuous operators) chosen as A →V_ℓℓ' (t) = V_ℓ (t) δ_ℓ, ℓ' for 
ℓ∈{0,1,2,…,N_bin-1}, representing the discrete fitness function,  
and B  → D_ℓℓ' = γ (δ_ℓ,ℓ'+1 - 2 δ_ℓ,ℓ' + δ_ℓ,ℓ'+1) as the discrete mutation matrix (note that the diagonal terms at ℓ =1 and ℓ=N_bin in the mutation matrix should be set to -γ). 
It is then straightforward to obtain the corresponding commutator matrix as 
C_ℓℓ' = γ (V_ℓ' - V_ℓ )(δ_ℓ,ℓ'+1+δ_ℓ,ℓ'-1). Following through the same steps leading to Eq. (<ref>), one can obtain the BCH approximation for the bnAb count in the discrete-bin model as 


    n_ℓ=0(t_f)  ≈      
        N_0 e^-λ t_f 
        e^∫_0^t_fdt V_ℓ=0(t)
    ×   e^γ∫_0^t_f dt (t_f-t) (V_ℓ=1(t) - 2V_ℓ=0(t) + V_ℓ=-1(t) )  
    ×   [e^t_f D_ℓ'ℓ” n_ℓ”(t=0)].

For the Gaussian fitness profile, optimizing this expression leads to the following transcendental equation

    exp1/2σ_V^2(t)/1/σ_V^2(t)-1 = 2γ(t_f-t)/1-2γ(t_f-t) .

This equation has a unique solution for 
2γ(t_f-t) ≥1/1+2e^-3/2 and no solution otherwise, defining a temporal window preceding the final time. It can be readily verified that the optimal width that solves Eq. (<ref>) decreases  over time until it reaches this final temporal window, beyond which the width remains fixed. Therefore, within the discrete-bin dynamics, the optimal fitness spread automatically reaches a finite minimum value, preventing the fitness amplitude from becoming singular, as discussed in Section <ref>.  



  
 

§ TIME EVOLUTION OF B CELL POPULATION UNDER DIFFERENT VACCINE PROTOCOLS
 
 



Here, we present further computational results on the time evolution of the B cell population in the one-dimensional shape space. 

Figure <ref> shows the bnAb count (top) and population peak's size (bottom) as functions of time. These results correspond to the same vaccination protocols and parameters as the left panel of Fig. <ref>. The graphs suggest that although positioning the vaccine-induced fitness peak optimally leads to the highest final bnAb count, the resulting population peak remains smaller compared to that achieved with the  χ^1,1/2_μ_0,σ_0(t) protocol. 
This observation is reinforced by the population snapshots in Fig. <ref>, indicating that while the x^*(t)  protocol effectively maximizes n_f(0), it may not fully optimize other aspects of the final B cell population, such as its total size or peak.  


In Fig. <ref>, we present the location of the population's peak as a function of time for the same vaccination protocols. Notably, only for the (χ^1/2,1_μ_0,σ_0,σ_min) protocol, the population peak appears to reach the target bin, whereas for the other two protocols, the peak remains away from x=0 at the final time.   
For the optimal vaccine center, and after an initial transient period, the population peak consistently trails x^*(t) with a constant lag. Conversely, with other choices, the population's distance from the vaccine center fluctuates over time. As discussed following Eq. (<ref>), such fluctuations in distance may either result in a slow movement of the population towards the target when the vaccine center remains too close to the population peak (blue) or lead to limited population growth when the vaccine center moves too far from the population (green). 

























 

