A Minimal Framework for Optimizing Vaccination Protocols
Targeting Highly Mutable Pathogens

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 [20, 17, 28]. Here we adopt a more coarse-grained approached based on the concept of shape-space representation [26, 27, 21]. In particular, BCRs and antigens are represented by points in a $d$-dimensional Euclidean space (Fig. 1). 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 $\bm{x}=(x_{1},x_{2},\ldots,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 [26, 29]. Following the approach outlined in Ref. [21], we characterize the binding free energy between a BCR and an antigen located in similarity bins $\bm{R}_{\mathrm{BCR}}$ and $\bm{R}_{\mathrm{Ag}}$ by their distance as

$$E_{\mathrm{bind}}\propto k_{\mathrm{B}}T\,||\bm{R}_{\mathrm{BCR}}-\bm{R}_{\mathrm{Ag}}||^{2}.$$

Consequently, the associated binding affinity is expressed as $(E_{\mathrm{thr}}-E_{\mathrm{bind}})$, with $E_{\mathrm{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\left[(E_{\mathrm{thr}}-E_{\mathrm{bind}})/{k_{\mathrm{B}}T}\right]$. 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 [30]. Furthermore, we assume that BCRs located closer to the origin ($\bm{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 $\bm{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 $\sigma_{\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.¹¹1We 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 [21]. 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. 1). 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(\bm{x})$, with its peak located at $\bm{x}_{V}$ (referred to as the vaccine center or focus) and its spread denoted as $\sigma_{V}$. It is worth noting that under the assumption of similar antigens, $\bm{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(\bm{x},t=0)\equiv n_{0}(\bm{x})$. As affinity maturation progresses, this distribution evolves over time to $n(\bm{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(\bm{x},t)$, across different bins in shape space, with corresponding time-dependent vaccine center and spread, $\bm{x}_{V}(t)$ and $\sigma_{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 $\gamma_{\bm{x},\bm{y}}$, corresponding to jumps between shape-space bins $\bm{x}$ and $\bm{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 $\lambda_{\bm{x}}$, representing their baseline apoptotic tendency within GCs when not receiving positive selection signals [5, 31].

The rules outlined above describe a generalized linear birth-death-mutation process and are encapsulated in the following Fock space master equation [32, 33]

	$\displaystyle\partial_{t}P(\{n\},t)=\quad$	$\displaystyle\sum_{\bm{x}}V_{\bm{x}}(t)\,\big{[}\mathbb{E}_{\bm{x}}^{-}-1\big{]}\big{(}n_{\bm{x}}P(\{n\},t)\big{)}$
	$\displaystyle\quad+$	$\displaystyle\sum_{\bm{x},\bm{y}}\gamma_{\bm{x},\bm{y}}\big{[}\mathbb{E}^{+}_{\bm{x}}\mathbb{E}_{\bm{y}}^{-}-1\big{]}\big{(}n_{\bm{x}}\,P(\{n\},t)\big{)}\,$
	$\displaystyle\quad+$	$\displaystyle\sum_{\bm{x}}\lambda_{\bm{x}}\,\big{[}\mathbb{E}^{+}_{\bm{x}}-1\big{]}\big{(}n_{\bm{x}}P(\{n\},t)\big{)}\,.$		(1)

Here, the BCR population is denoted by $\{n\}=\{n_{1},n_{2}\,\ldots\}$ where $n_{i}$ is a non-negative integer representing the occupation number of the $i$th bin, and $P(\{n\},t)$ is the probability of having that specific population at time $t$. On the right-hand side of Eq. (II.1), we introduce the operators $\mathbb{E}_{\bm{x}}^{\pm}$, defined as

$$\mathbb{E}_{\bm{x}}^{\pm}\big{(}n_{\bm{y}}P(\{n\},t)\big{)}\equiv(n_{\bm{y}}\pm\delta_{\bm{x},\bm{y}})\,P(\{n\}\pm\mathbb{I}_{\bm{x}})\,,$$

where $\bm{x}$ and $\bm{y}$ are bin indexes, $\delta_{\bm{x},\bm{y}}$ is the Kronecker delta, and $\{n\}\pm\mathbb{I}_{\bm{x}}$ denotes a population that differs from $\{n\}$ by one more or less B cell residing in bin $\bm{x}$.

Each line on the right-hand side of Eq. (II.1) 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. (II.1) would provide a probabilistic description of the entire BCR population in the high-dimensional shape space.

The abstract model described above offers a simplified perspective on affinity maturation. However, deriving insights from the resulting master equation (Eq. (II.1)) 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 $\expectationvalue{n_{\bm{x}}(t)}=\sum_{\{n\}}n_{\bm{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 A for derivation)

$\displaystyle\partial_{t}n(\bm{x},t)=D\,\nabla^{2}n(\bm{x},t)+\left[V(\bm{x},t)-\lambda\right]\,n(\bm{x},t)\,,$

(2)

where, for notational simplicity, we retain the symbol $n$ to represent the average population density in the continuum limit of $\bm{x}$. We also assume the death rate $\lambda$ and the diffusion coefficient $D$ (obtained from the underlying mutation rates via Eq. (31)) are constants independent of $\bm{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 $\bm{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. [18], where antigens are different enough that bnAb sequences are always the most fit.

Solving the mean-field dynamics (Eq. (2)) requires specifying the initial germline distribution, which typically has little overlap with the target bin in shape space. We thus assume $n_{0}(\bm{x})$ is centered at a distant location $\bm{\mu}_{0}$ and has a width $\sigma_{0}$ that is an indication of the diversity of the germline BCRs (see Fig. 1). In the spirit of the central limit theorem, we model $n_{0}(\bm{x})$ by a normal distribution

$\displaystyle n_{0}(\bm{x})=N_{0}\,\dfrac{e^{-\frac{(\bm{x}-\bm{\mu}_{0})^{2}}{2\sigma_{0}^{2}}}}{\left(2\pi\sigma_{0}^{2}\right)^{d/2}}\equiv N_{0}\mathcal{G}(\bm{x}-\bm{\mu}_{0},\sigma_{0}^{2})\,,$

(3)

where $\mathcal{G}$ 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 $\{\bm{x}_{V}(t),\sigma_{V}(t)\}$, that maximizes the final bnAb population denoted by $n(\bm{x}=0,t_{f})\equiv n_{f}(0)$, while adhering to a constrained total vaccine dose, i.e.

$$\int_{0}^{t_{f}}dt\int d^{d}\bm{x}\,V(\bm{x},t)=\text{constant}\,.$$

(4)

Equation (2) resembles the imaginary-time Schrodinger equation with a time-dependent potential [34]. 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:

1. 1.

   We interpret Eq. (4) 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. (4)) and the amount of the administered antigen. However, $V(\bm{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. [5, 35, 36, 37]). 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. (II.1).

2. 2.

   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. (2), along with the vanishing drift term (see Appendix A). 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 [36], potentially leading to long-ranged steps in the coarse-grained mean-field dynamics.

We begin constructing the path integral formulation by rearranging Eq. (2) for an infinitesimal timestep $\delta t$:

	$\displaystyle n(\bm{x},t+\delta t)$	$\displaystyle=\{1+\delta t\,D\nabla^{2}+\delta t[V(\bm{x},t)-\lambda]\}n(\bm{x},t)$
		$\displaystyle=e^{\delta tD\nabla^{2}}e^{\delta t[V(\bm{x},t)-\lambda]}n(\bm{x},t)+O\left(\delta t^{2}\right)\,.$		(5)

Equation (III.1) provides a formal method for evolving the population profile in infinitesimal timesteps. By repeatedly applying Eq. (III.1) to evolve the population over time, the final population $n(\bm{x},t_{f})\equiv n_{f}(\bm{x})$ can be expressed by the formal time-ordered product [34]

$\displaystyle n_{f}(\bm{x})=\lim_{M\to\infty}\prod_{j=0}^{j=M}e^{\delta tD\nabla^{2}}\,e^{\delta t[V(\bm{x},j\delta t)-\lambda]}\,n_{0}(\bm{x})\,.$

(6)

Here, we have defined $\delta t=\frac{t_{f}}{M}$ with $M$ as the number of time increments, and $V_{j}(\bm{x})$ denotes the vaccine-induced fitness function at time $j\delta t$. The products in Eq. (6) act in the operator sense:²²2This product formula is also known as the first-order Suzuki–Trotter decomposition [38, 39] and it forms the basis for the path integral formulation of quantum mechanics [34], quantum simulators [40, 41], and the so-called operator splitting methods [42]. starting from the germline population $n_{0}$, at the $j$th infinitesimal timestep $\delta t$, the B cell population is locally expanded/contracted according to the vaccine operator $e^{\delta t[V(\bm{x},j\delta t)-\lambda]}$, followed by diffusion broadening under the operator $e^{\delta tD\nabla^{2}}$.

The effect of the diffusion operator on a test function $\eta(\bm{x})$ is given by the convolution

$$e^{\delta tD\nabla^{2}}\eta(\bm{x})=\int d^{d}\bm{x}^{\prime}\mathcal{G}(\bm{x}-\bm{x}^{\prime},2D\delta t)\,\eta(\bm{x}^{\prime}).$$

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. (6) and integrating over all intermediate points, the final population can be more explicitly expressed as

$$\int\prod_{j=0}^{M}d^{d}\bm{x}_{j}\,\mathcal{G}(\bm{x}_{j+1}-\bm{x}_{j},2D\delta t)\,e^{\delta t(V(\bm{x}_{j},j\delta t)-\lambda)}\,n_{0}(\bm{x}_{0})\,,$$

with $\bm{x}_{M+1}=\bm{x}$. This expression is then rearranged into the standard path integral form [34, 43]

$\displaystyle n_{f}(\bm{x})=e^{-\lambda t_{f}}\int d^{d}\bm{x}_{0}\,n_{0}(\bm{x}_{0})\int_{\bm{r}(0)=\bm{x}_{0}}^{\bm{r}(t_{f})=\bm{x}}\mathcal{D}\bm{r}(u)\,e^{-S[\bm{r},\dot{\bm{r}}]}\,,$

(7)

where $\mathcal{D}\bm{r}\propto d\bm{x}_{1}\ldots d\bm{x}_{M}$ stands for the measure of all paths connecting the endpoints, each weighted according to its ‘action’

$\displaystyle S[\bm{r},\dot{\bm{r}}]=\int_{0}^{t_{f}}du\,\left(\frac{\dot{\bm{r}}^{\,2}(u)}{4D}-V(\bm{r}(u),u)\right)\,.$

(8)

Note that Eq. (7) indicates that evolutionary paths passing through regions with high vaccine fitness tend to have larger weights in the path integral expression.

Equation (7) 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 (3) into Eq. (7) and then apply the linear coordinate transformation $\bm{x}\to\bm{x}+\bm{x}^{*}(t)$, where

$\displaystyle\bm{x}^{*}(t)\equiv\frac{\bm{\mu}_{0}}{1+\frac{\sigma_{0}^{2}}{2Dt_{f}}}\left(1-\frac{t}{t_{f}}\right)\,.$

(9)

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

$\displaystyle\dfrac{n_{f}(0)\big{\rvert}_{\bm{\mu}_{0}=\bm{\ell};V}}{n_{f}(0)\big{\rvert}_{\bm{\mu}_{0}=0;\widetilde{V}}}=\exp{-\frac{\bm{\ell}^{2}}{2(\sigma_{0}^{2}+2Dt_{f})}}\,.$

(10)

In this equation, the numerator on the left-hand side denotes the final population at $\bm{x}=0$ (i.e., the number of final bnAbs) given that the normally-distributed initial population (3) is centered at $\bm{\mu}_{0}=\bm{\ell}$ and evolves under the vaccine fitness $V(\bm{x},t)$. The denominator, on the other hand, represents the final bnAb count for a germline population that is centered at $\bm{\mu}_{0}=0$ and is evolved under the transformed fitness $\widetilde{V}$ defined as

$\displaystyle\widetilde{V}(\bm{x},t)\equiv V\left(\bm{x}+\bm{x}^{*}(t),t\right)\,.$

(11)

Therefore, the exact relationship (10) 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 $\widetilde{V}$, assuming the germline is already centered at the target (see top panel of Fig. 2).

The crucial observation here is that the right-hand side of Eq. (10) 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 $\{\bm{x}_{V}(t),\sigma_{V}(t)\}$), must coincide with the maximum of the denominator as $\widetilde{V}$ is varied. As a result, Eq. (10) effectively converts the maximization of $n_{f}(0)\big{\rvert}_{\bm{\mu}_{0}=\bm{\ell};V}$ into the maximization of $n_{f}(0)\big{\rvert}_{\bm{\mu}_{0}=0;\widetilde{V}}$. This mapping simplifies the optimization problem: since for $\bm{\mu}_{0}=0$, the starting germline population is already peaked at, and symmetric around, the target bin $\bm{x}=0$, it makes sense that for $\widetilde{V}$ to be optimal, it should also be peaked at the target and symmetric with respect to it at all times (see Fig. 2). Using this argument, and recalling that the maxima of the numerator and denominator in Eq. (10) coincide, we then conclude that

$$\bm{x}_{V}(t)\big{\rvert}_{\text{optimal}}=\bm{x}^{*}(t)=\bm{\mu}_{0}~{}\frac{2D(t_{f}-t)}{{2Dt_{f}}+{\sigma_{0}^{2}}}\,.$$

(12)

It is worth noting that for a narrow germline population ($\sigma_{0}\ll\sqrt{Dt_{f}}$), the optimal vaccine center given by Eq. (9) simply starts close to the naive population at $\bm{\mu}_{0}$ and reaches the target state ($\bm{x}=0$) at the final time. However, for a wider initial population ($\sigma_{0}\gtrsim\sqrt{Dt_{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. (7). Intuitively, one needs to position the vaccine fitness over time such that it gradually moves the germline population from $\bm{x}=\bm{\mu}_{0}$ towards the target state at $\bm{x}=0$, while also remaining sufficiently close to the peak of the existing B cell population to avoid extinction by apoptosis. If $V(\bm{x},t)$ is positioned too close to the population $n(\bm{x},t)$, it may not optimally move the population to reach $\bm{x}=0$ within the finite time interval $t_{f}$; conversely, if $V$ is centered too far from the peak of $n(\bm{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.  (9)) aligns with this intuitive understanding and also demonstrates how the germline width $\sigma_{0}$ influences the optimal position of the vaccine fitness. In addition, Eq. (12) 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 IV and Appendix D.

As discussed above, Eq. (10) maps the problem of maximizing bnAbs to finding the optimal modified fitness function $\widetilde{V}$ that maximizes $n_{f}(0)\big{\rvert}_{\bm{\mu}_{0}=0;\widetilde{V}}$. Upon replacing $V$ by $\widetilde{V}$, Eq. (2) describes the expansion of the transformed germline population centered at $\bm{x}=0$ under the modified fitness $\widetilde{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. (6), we rearrange the operators in the time-ordered product formula using a second-order truncation of the celebrated Baker-Campbell-Hausdorff (BCH) formula [39], namely

$$e^{\epsilon A}e^{\epsilon B}=e^{\epsilon B}e^{\epsilon A}e^{\epsilon^{2}[A,B]}\left(1+O(\epsilon^{3})\right),$$

(13)

where $A$ and $B$ are generally non-commuting operators, and $C\equiv[A,B]=AB-BA$. In this formula, we assign $\epsilon\to\delta t$, $A\to D\nabla^{2}$, and $B\to\widetilde{V}_{j}(\bm{x})$, and then utilize this formula to reorder the alternating operator product in Eq. (6) such that all vaccine fitness operators $e^{\delta t\,\widetilde{V}_{j}}$ are grouped on the left in the product, all the diffusion operators $e^{\delta tD\nabla^{2}}$ are collected on the right, and the commutator terms that result from swapping these operators appear in between. The apoptosis operator $e^{-\lambda t_{f}}$ commutes with all the other operators.

To achieve such an ordering of the operators, we start by moving $e^{\delta t\widetilde{V}_{0}}$ through the $M$ diffusion operators $e^{\delta tD\nabla^{2}}$ on its left in Eq. (6), resulting in a net commutator term $e^{(\delta t)^{2}MC_{0}}$, where the commutator operator reads

$$C_{0}=[D\nabla^{2},\widetilde{V}_{0}(\bm{x})]=2D\widetilde{V}_{0}^{\prime}(\bm{x})\nabla+D\widetilde{V}_{0}^{\prime\prime}(\bm{x})\,.$$

(14)

Subsequently, the next fitness operator, $e^{\delta t\widetilde{V}_{1}}$, needs to pass through $(M-1)$ diffusion operators on its left, thereby creating a commutator contribution $e^{(\delta t)^{2}(M-1)C_{1}}$. This process continues for all fitness operators, and in the limit of small increments $\delta t=\frac{t_{f}}{M}\to 0$, we obtain³³3Each truncation of the BCH formula introduces an error of order $(\delta t)^{3}$. It can be shown that the number of such error terms grows at least as $\sim 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.

	$\displaystyle n_{f}\big{\rvert}_{\bm{\mu}_{0}=0;\widetilde{V}}\approx$	$\displaystyle N_{0}\quad e^{-\lambda t_{f}}e^{\int_{0}^{t_{f}}dt\,\widetilde{V}(\bm{x},t)}$
		$\displaystyle\times e^{D\int_{0}^{t_{f}}dt\,(t_{f}-t)\,\left(2\widetilde{V}^{\prime}(\bm{x},t)\nabla+\widetilde{V}^{\prime\prime}(\bm{x},t)\right)}$
		$\displaystyle\times e^{t_{f}D\nabla^{2}}\,\mathcal{G}(\bm{x}_{0},\sigma_{0}^{2})\,.$		(15)

Unlike the exact expression in Eq. (6), where the fitness and mutation (diffusion) alternate in the operator product, Eq. (III.2) suggests that $n_{f}$ can be approximated by first diffusing the starting profile $n_{0}$ under the (full) diffusion operator $e^{t_{f}D\nabla^{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\nabla^{2}}$ on the Gaussian initial population widens its profile and yields $\mathcal{G}(\bm{x}_{0},\sigma_{0}^{2}+2Dt_{f})$, which is independent of the fitness profile. The symmetry of this updated profile with respect to $\bm{x}_{0}=0$ still implies a similar symmetry for the optimal $\widetilde{V}(\bm{x},t)$, resulting in $\widetilde{V}^{\prime}(\bm{x}=0,t)=0$. Consequently, Eq. (III.2) simplifies to

	$\displaystyle n_{f}(0)\big{\rvert}_{\bm{\mu}_{0}=0;\widetilde{V}}\approx\quad$	$\displaystyle N_{0}e^{-\lambda t_{f}}e^{\int_{0}^{t_{f}}dt\widetilde{V}(0,t)}$		(16)
	$\displaystyle\times$	$\displaystyle e^{D\int_{0}^{t_{f}}dt(t_{f}-t)\widetilde{V}^{\prime\prime}(0,t)}\,\mathcal{G}(0,\sigma_{0}^{2}+2Dt_{f})\,.$

For any choice of the vaccine-induced fitness profile, Eq. (16) 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

$\displaystyle V_{\mathcal{G}}(\bm{x},t)\equiv A_{V}\,\mathcal{G}\left(\bm{x}-\bm{x}_{V}(t),\sigma_{V}(t)^{2}\right)\,.$

(17)

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. (17) also satisfies the total-dose constraint given by Eq. (4). However, it is important to note that Eq. (17) 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 [44, 45, 46]. We defer the investigation of such scenarios to future studies.

By selecting the optimal vaccine center according to Eq. (12), the Gaussian fitness of Eq. (17) transforms according to Eq. (11), resulting in the modified fitness $\widetilde{V}_{\mathcal{G}}(\bm{x},t)=A_{V}\mathcal{G}(\bm{x},\sigma_{V}^{2}(t))$, which is still a Gaussian profile but now centered at the origin. Substituting $\widetilde{V}_{\mathcal{G}}$ into Eq. (16), we obtain

	$\displaystyle n_{f}(0)\big{\rvert}_{\bm{\mu}_{0}=0;\widetilde{V}_{\mathcal{G}}}\approx\frac{N_{0}e^{-\lambda t_{f}}}{\sqrt{2\pi(\sigma_{0}^{2}+2Dt_{f})}}$		(18)
	$\displaystyle\qquad\times\exp{\int_{0}^{t_{f}}dt\,\left(\frac{A_{V}}{\sqrt{2\pi}\sigma_{V}(t)}\left(1-\frac{D(t_{f}-t)}{\sigma_{V}^{2}(t)}\right)\right)}\,.$

To maximize this expression, we set its functional derivative with respect to $\sigma_{V}$ to zero. Assuming $\sigma_{V}(t)\geq\sigma_{\min}$ at all times (see the discussion below), we finally obtain

$\displaystyle\sigma_{\mathrm{BCH}}(t)=\mathrm{max}\left\{\sqrt{3D(t_{f}-t)}\,,\sigma_{\min}\right\}$

(19)

In deriving Eq. (19), we have assumed that $\sigma_{V}(t)\geq\sigma_{\min}$, which prevents $V_{\mathcal{G}}$ from reaching arbitrarily large values by imposing a maximum fitness value $A_{V}/(\sqrt{2\pi}\sigma_{\min})$. As discussed in Section II, 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 C, 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 $\sigma_{\min}$ encompasses all of these factors.

Equation (19) can be understood by noting that within a temporal window of duration $\Delta t$, diffusion affects a range of order $\sqrt{D\Delta 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 $\widetilde{V}$ should at most cover BCRs within a range of order $\sqrt{D(t_{f}-t)}$ while remaining peaked at $\bm{x}=0$. Once the vaccine spread reaches the smallest possible value $\sigma_{\min}$, it should remain at that value until the final time. In a qualitative sense, Eq. (19) 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. (19), 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 (16) depends only on time integrals of $\widetilde{V}(0,t)$ and $\widetilde{V}^{\prime\prime}(0,t)$, making it a strictly local function of the transformed vaccine fitness $\widetilde{V}$ at $\bm{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 $\bm{x}$. Since the mean-field Eq. (2) 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 [11, 13], 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 II.2, our approach differs from Refs. [18, 17] 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. (19) bears resemblance to the findings of Ref. [18]– where the optimal fitness function appeared narrower in subsequent vaccinations–we approach the interpretation of Eq. (19) with caution. In the next section, we will compare the performance of protocols employing a width $\sigma_{\mathrm{BCH}}$ with those utilizing a narrow fitness width $\sigma_{\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.