A predictive systems vaccinology framework enables rational optimization of MVA-based vaccines

Reviewer #1 (Public review):

Summary:

The manuscript "A predictive systems vaccinology framework enables rational optimization of MVA-based vaccines" by Deman and co-workers presents an approach to use Boolean models for the optimization of MVA for vaccinations. Different Boolean models are derived/inferred to perform in silico testing, e.g., of knock-outs.

Strengths:

The optimization of vaccine platforms is very important, and model-based approaches have proved a powerful framework for in silico testing. As far as I'm aware, this is the first time a comprehensive Boolean model is used for this. The authors make an effort to inform this model from available information and experimental data, using state-of-the-art calibration pipelines.

Weaknesses:

(1) Lines 154-158: "Because certain biological processes represented in KEGG (e.g., phosphorylation or ubiquitination) do not have direct logical equivalents, this conversion of signaling pathways into a Boolean network can lead to information loss and disconnection of nodes from the rest of the network. To mitigate this issue, we reconnected isolated nodes back to the main structure using oriented protein-protein interaction (PPI) data from 69, thereby restoring connectivity while preserving directionality of regulation." It is not clear to me how the reconnection addresses the described issue that not all processes can be represented in the selected modelling framework. In this context, I would also appreciate it if the authors could clarify the meaning of your states. Is it the presence of a protein (relating to low/high abundance), the activation status (relating to low/high phosphorylation), or a combination? Depending on this, different Boolean representations should be chosen, and different process information can be used.

(2) Lines159-160: "To enhance immediate readability and interpretability, we connected the resulting network with the corresponding cellular population abundances analyzed by cytometry in the samples." I would appreciate it if the authors could clarify how the cellular layer and the population layers were connected. Is this related to proliferative potential?

(3) Line 168++: It is unclear to me which parts of the Boolean network described in the section "Boolean naïve network construction" have been calibrated. Among other things, it would also be interesting to know how many logical expressions were changed by ZhegAlCal compared to the naive model and how these expressions were selected. Is there a regularization aiming to minimize the number of changes? In this context, I would also appreciate a clarification of the data processing. The current text mentions a 20% change compared to baseline, a threshold of 0.05, and a 2-means clustering strategy, yet it is unclear how they interact to obtain the binarized training and validation data.

(4) Line 267++: The model constructed by the authors describes cell-level processes in infected cells. Yet, the data used in the study - which have previously been published in reference 47 - seem to rather capture population averages over heterogeneous, partially non-infected cells. It is unclear to me why / how this can be compared. I would appreciate a clarification, potentially including a more detailed description of the employed datasets.

(5) Lines 758-759: "The networks generated and analysed during this study are publicly available in the CellCellective repository (MVA 3 pathways, MVA 6 pathways, YF17D)." I searched for the research but did not find it. In my opinion, it would be important to make the models as well as the implementations for calibration, etc. available. Without this, value and reproducibility are limited. I would encourage the authors to provide a detailed human-readable model description in the supplement.

(6) Lines 783-785: The GO analysis seems to be performed in comparison to the human genome. Yet, the model contains only 200 nodes, so a substantially reduced fraction. I was wondering if this was considered in the analysis process and if the authors checked how often the enrichments for multiple pathways were driven by the same genes.

(7) Figure 3: It appears as if the number of considered "network updates" was set to 10 (0 to 9) and that this somehow maps to the experimental time. Yet, the experimental observation times are far from uniform.

Reviewer #2 (Public review):

Summary

Boolean network modeling is more commonly used in cancer and developmental biology than in vaccine research. Deman et al. apply this framework to a practical vaccinology problem: they aimed to build a mechanistic, executable computational framework capable of both explaining and predicting how the early innate immune response to the MVA vaccine changes when specific viral genes are altered, with the longer-term goal of using that framework to guide the rational design of improved MVA-based vaccines. The authors aimed to: (i) construct and calibrate a Boolean network model of the MVA-induced immune response against real longitudinal non-human primate (NHP) data; (ii) test whether the calibrated model, without being fit to this new data, could reproduce the outcomes of several previously published MVA gene-deletion mutants; and (iii) compare this MVA model to an analogous model of the well-established YF-17D yellow fever vaccine, in the hope of identifying specific molecular targets that could reorient the MVA response toward the durable, single-dose protection YF-17D is known to provide.

In pursuit of these aims, the authors construct an executable Boolean network of the innate immune response to the MVA vaccine by merging three KEGG pathways (cytosolic DNA-sensing, apoptosis, NF-κB signaling) with cell-population data, and calibrate it against a small NHP dataset (n=3 macaques, 7 time points; Rosenbaum et al., ref. 47). They show the calibrated network reproduces 86-87% of the observed binarized states, and that forcing the network to mimic known MVA deletion mutants (e.g., the triple mutant ΔC6L/ΔK7R/ΔA46R) reproduces qualitative features (e.g., early IFN-β, TNF-α, and IL-6 upregulation) reported in independent published mouse and cell-line studies. They then build a second, six-pathway "consensus" network shared between MVA and YF-17D, compare the two networks' topology and dynamics, and use this comparison, together with a graph-theoretic search for "highly effective" signaling paths, to propose two previously untested MVA deletion mutants (ΔK7R and ΔF17R) predicted to shift the MVA response toward more YF-17D-like features.

Strengths

The overall workflow (Figure 1) is clearly described. The calibration approach-binarizing longitudinal cellular and transcriptomic data and fitting network trajectories with the ZhegAlCal algorithm (a Zhegalkin-polynomial/SAT-solving-based method for fitting Boolean trajectories to binarized time-series data)-appears to be a defensible, well-reasoned way to translate a literature-derived network into real kinetic data.

The retrospective validation against independent published MVA mutants (deletions in C6L, K7R, A46R, and N2L, and separately an A21L point-mutant with three alanine substitutions rather than a deletion) is a genuine strength: the model's qualitative behavior (upregulation of IFN-β, TNF-α, IL-6; limited change in RIG-I) aligns with what those studies reported. We also appreciated that the authors report instances of partial disagreement alongside their successes (e.g., CCL5/RANTES) - that kind of candor about where the model doesn't quite line up is exactly what gives the parts that do line up more credibility.

The static topological analysis - hub identification, "determinative power" and "vertex betweenness" (two complementary measures of how much a node's state constrains, or lies on paths between, the rest of the network), and "effective graphs" (a measure of how deterministically an edge's regulator sets its target's state) - is a thoughtful use of graph theory to complement the dynamic simulations.

The authors also clearly discuss the Boolean formalism's core approximation, i.e., that binary on/off states can dilute real but subtle quantitative differences (lines 679-688), and that this work was based on modeling blood-only responses rather than those responses that occur at the vaccination site or draining lymph nodes (lines 709-714).

Weaknesses

A few things gave us pause as we read, which we raise here in the spirit of strengthening what already strikes us as a promising framework.

The MVA/YF-17D comparison starts from two vaccines that are already known to differ substantially. The manuscript uses the divergence between the MVA and YF-17D Boolean networks as an entry point for identifying "MVA optimization" opportunities, but MVA and YF-17D are, on their face, very different vaccine platforms. MVA is a non/limited-replicating DNA poxvirus vector, sensed mainly through cytosolic DNA/cGAS-STING pathways, dosed intradermally, and typically requiring two doses for optimal protection. YF-17D, by contrast, is a live, replicating, attenuated RNA flavivirus, sensed through multiple TLR/RIG-I pathways, and given as a single subcutaneous dose that confers durable, often lifelong, protection (see refs 21, 31, 38-44 in the manuscript). Given this, it's not surprising that the authors themselves report "almost opposite behaviors" for several core cell populations - classical monocytes, B cells, NK cells, and CD4/CD8 T cells - between the two calibrated networks (lines 647-657).

This stated motivation made us question how much of that divergence reflects a real, actionable difference in vaccine-induced immune programming (the paper's implicit premise) versus differences in virus biology, dosing route, or study/technical design (different sampling schedules, microarray vs. RNA-seq, n=3 vs. n=12 animals) that a Boolean network comparison can't easily tease apart. To their credit, the authors' Discussion is candid on this point-stating directly that "no clear optimization strategies for the MVA viral vector came to mind from the comparison" (lines 664-666)-a useful signal that the comparison's direct yield was limited. The two candidate mutants that emerged instead came from intersecting the model's high-impact nodes with a pre-existing, literature-curated list of MVA immunomodulatory genes (yielding five candidate deletions), which were then individually simulated and narrowed down to the two, ΔK7R and ΔF17R, that produced notable changes - not from the MVA/YF-17D comparison alone.

The motivation for choosing Boolean modeling over ODE/PDE approaches could be clearer. The choice is motivated mainly by precedent - the authors note it "has rarely been used to model vaccine-induced immune responses, in contrast to statistical modeling or ordinary differential equation (ODE)-based modeling" (lines 100-105) - and by practical considerations, such as not requiring kinetic parameters and being tractable at the scale of networks with hundreds of nodes. What we found ourselves wanting was a more explicit account of the trade-off: ODE and PDE models already simplify the true, spatially resolved, continuously varying underlying biology, and Boolean modeling is a further simplification on top of that. A clearer statement of why this additional simplification is acceptable, or even preferable, for this application (for example: the scale of the curated network, the absence of measured rate constants for most edges, or the interpretability of discrete states) would greatly improve this work.

We found the manuscript dense and long relative to the size of its central, generalizable findings. This is a presentation issue rather than an evidentiary one. The Results section narrates GO-enrichment interpretation update-by-update for four separate network trajectories (unperturbed MVA, perturbed MVA mutants, the MVA arm of the consensus network, and YF-17D), much of which restates what the (extensive) supplementary figures already show.

The paper's only prospective predictions are experimentally untested. The two new mutants proposed at the end of the paper (ΔK7R, ΔF17R) are in silico predictions only; they haven't been constructed or tested experimentally in this study. We read the title's claim of a "predictive" framework and the paper's "rational design" framing as best describing a hypothesis-generation tool validated by retrodiction of previously published phenotypes, rather than a demonstration that these two newly proposed mutants will behave as predicted in vivo.

Did the authors achieve their aims, and do the results support their conclusions?

Looking at each aim in turn, the picture that emerges is mixed. The first aim - building and calibrating a Boolean network model of the MVA-induced immune response - is convincingly achieved: the calibrated network fits the underlying data well (86-98% of binarized states correctly reproduced, depending on which of the two networks is considered), and its successive states correspond to biologically sensible processes (early chemotaxis, then T-cell activation, then antiviral/ROS-related signatures) that track what is independently known about the innate response to poxvirus vaccination. The second aim-showing the calibrated model can reproduce, without being fit to it, the outcomes of previously published MVA mutants-is also substantially achieved, with the caveat noted above that agreement is qualitative and directional rather than exact, and not uniform across every marker tested.

The third aim is where the results, in our reading, support the paper's conclusions least well. The stated purpose of comparing the MVA and YF-17D networks was to identify actionable strategies for reorienting MVA's response, but the authors themselves report that the comparison alone yielded no clear optimization strategy (lines 664-666); as noted above, the two candidates that are ultimately proposed came instead from that separate gene-list intersection and simulation step - not from the YF-17D comparison in the more direct way the framing implies. Given that, the paper's headline conclusion - that this framework "enables rational optimization of MVA-based vaccines" - reads to us as only partially supported by what is actually shown: the framework is well demonstrated as a tool for capturing and reproducing known immune biology, but its capacity to prospectively guide the design of a better vaccine remains, at this point, an untested hypothesis rather than a demonstrated result.

Likely impact and utility to the community:

The most durable contribution of this paper, regardless of how the two specific candidate mutants eventually fare in the laboratory, strikes us as methodological: the explicit workflow for merging curated signaling pathways into a large executable Boolean network and calibrating it against longitudinal experimental data (building on the authors' own previously published calibration method, reference 75) is clearly described and, together with the deposition of the calibrated networks on the public CellCollective platform, should be usable by other groups working on other vaccines or viral vectors. That reusability is a genuine and useful contribution to the systems-vaccinology toolkit, independent of whether MVA specifically benefits from it.

Its more immediate, practical utility is harder to gauge from the paper alone. For vaccine developers specifically interested in MVA, the value of this work currently lies in the two testable hypotheses it generates (ΔK7R, ΔF17R) rather than in validated design guidance, since neither candidate has been built or tested here. It's also worth flagging that the framework's generalizability beyond MVA and poxviruses is untested within this paper - the approach is demonstrated for one vector and one comparator vaccine, so readers working on other vaccine platforms may want to treat it as a promising template to adapt and validate for their own systems, rather than as a result that has already been shown to transfer.