Introduction

The poxvirus Modified Vaccinia virus Ankara (MVA) is a highly attenuated strain of vaccinia virus (VACV) that does not fully replicate in mammalian cells, and is widely used as a safe and versatile viral vaccine vector. Originally developed as a third-generation smallpox vaccine, MVA is now licensed for use against mpox to face the 2022 outbreak ^1–3 and as a recombinant booster vaccine in Ebola vaccination regimens ⁴, and serves as a versatile vector for developing vaccines against numerous infectious diseases and cancers ^5–7. MVA vaccines are known for their strong immunogenicity and protective efficacy against poxviruses, as they elicit robust B- and T-cell responses while maintaining an excellent safety profile ^8–11. The induction of long-lasting antigen-specific adaptive immunity by vaccination is known to be strongly orchestrated and modulated by the early innate immune responses ^12–14 For MVA, specifically, early defense mechanisms have been shown to strongly affect long-term antibody and T-cell responses in preclinical models ^15,16 and in clinical trials in human volunteers ¹⁷. Nevertheless, its immunogenicity should be improved, as a two-dose regimen is necessary to induce comparable protection against mpox in macaques to a single-dose strategy with a live VACV vaccine ¹⁸, and the neutralizing anti-mpox antibody response wanes rapidly, particularly in smallpox vaccine-naïve individuals ^19,20. Furthermore, the molecular and cellular mechanisms linking early innate sensing of MVA to durable adaptive immunity remain only partially understood, limiting efforts to rationally optimize this vaccine vector.

The MVA genome, sequenced in 1998 ²¹, lacks 15% of the parental VACV genome, including several viral immune evasion and host-range genes. However, multiple viral immunomodulatory genes still remain conserved in the MVA genome and actively modulate and often inhibit host innate pathways ²². Over the past two decades, targeted deletions of such genes have yielded MVA recombinants with enhanced immunogenicity ^23–29, demonstrating that MVA vector engineering can reshape host immune responses. Yet, identifying which viral or host pathways should be modulated to achieve optimal immunogenicity remains largely empirical, highlighting the need for predictive frameworks capable of linking molecular mechanisms to system-wide immune outcomes.

The advent of systems vaccinology, an interdisciplinary approach that leverages computational and experimental methods to understand and predict vaccine-induced immune responses, opens the door to a more rational definition of vaccine optimization strategies ³⁰. To guide the rational search for optimized MVA-based vaccines, one promising approach is to mimic the characteristics of a highly effective vaccine-induced response. In this context, the live-attenuated YF-17D yellow fever vaccine represents a compelling benchmark. Considered as a paradigm of vaccine success, YF-17D combines an excellent safety profile (despite rare, around 1/250000 administrations, but sometimes severe adverse events) with lifelong protection after a single dose ³¹. Moreover, it has been emphasized the learning potential from YF-17D, especially as its mechanisms of action have been progressively elucidated through systems biology approaches ³². To date, such studies have provided detailed descriptions of the YF-17D-induced response and identified mechanisms of action and predictive markers of immunogenicity and protective efficacy ^33–36, as well as unique features that distinguish it from other vaccines ³⁷. YF-17D infects dendritic cells (DCs) and stimulates multiple Toll-like receptor (TLR) pathways ^38–44. Building on this knowledge, nanoparticles containing antigens combined with ligands for one or more TLRs were designed and evaluated in mice. The simultaneous stimulation of multiple TLR pathways was shown to induce stronger adaptive responses than single-pathway activation ⁴⁵. Similarly, the immunogenicity of a dengue vaccine candidate was enhanced by adjuvating the envelope antigen with the combination of TLR agonists mimicking the best immune signatures of early responses to YF-17D ⁴⁶. These examples illustrate how mechanistic understanding ofYF-17D can inform rational vaccine design.

In ⁴⁷, we conducted a non-human primate (NHP) study to characterize early immune responses following MVA vaccination. Using cytometry, transcriptomics, and cytokine profiling combined with systems biology approaches, including pathway enrichment, dimensionality reduction, and co-expression network analyses, we delineated the cascade of cellular and molecular events initiated at the vaccine injection site, in draining lymph nodes, and in blood. We linked the magnitude and quality of early innate effector responses with the orientation of subsequent adaptive immunity and highlighted a potential role for myeloid-derived suppressor cells (MDSCs). Building on these findings, dynamic mechanistic modeling represents a logical next step to extend our knowledge and understanding and improve our capacity to optimize MVA-induced immune responses. By enabling the simulation of MVA- and YF-17D-elicited immune dynamics over time, such models could provide insight into shared and distinct mechanisms and thereby guide the rational optimization of MVA vaccination strategies.

Boolean modeling is a dynamic and mechanistic formalism that translates the continuous expressions and dynamics of the underlying biological systems into binary states and qualitative changes ⁴⁸. It has been successfully applied to diverse contexts, ranging from cellular expression and plasticity ^49–51 to the pathogenesis of cancer and infectious diseases and potential therapeutic interventions ^52–56, but it has rarely been used to model vaccine-induced immune responses, in contrast to statistic modelling or ordinary differential equation (ODE)-based modelling ^57–64 Because Boolean formalism implies a substantial approximation with regard to the underlying biological system, it requires careful validation against experimental observations ^49,53,65. Importantly, to avoid overfitting and actually to challenge computational models, the validation should be performed in conditions outside the scope of the data used to build the models ⁶⁶.

In this work, we applied the Boolean formalism to the modeling of vaccine-induced immune responses. We first build a Boolean network of the MVA-induced immune response by integrating literature-curated existing knowledge with context-specific experimental data from our MVA preclinical study in NHPs ⁴⁷. The resulting model precisely recapitulated the MVA-induced immune response with high fidelity. To validate the underlying mechanisms and assess its extrapolation capacity, we leveraged published immune responses data from genetically modified MVA vectors ^23,25,28. By comparing perturbed simulation outcomes with their independent experimental observations, we demonstrated strong concordance, supporting the validity of both the model and, by extension, the modeling strategy, as well as the relevance and usefulness of Boolean modeling in vaccinology. We then exploited the executable nature of Boolean networks in both supervised and unsupervised modes to discover strategies for enhancing MVA immunogenicity. We first evaluated in silica two novel MVA mutants derived from expert knowledge by perturbing the model of response to MVA. Next, we built two comparable Boolean models of MVA- and YF-17D-induced immune responses in NHPs and performed static and dynamic comparative analyses on these networks to identify shared and divergent regulatory features. This comparative framework enabled the virtual identification and evaluation of molecular targets to reorient MVA-induced immune responses toward YF-17D-like immune signatures.

Materials and methods

The Boolean modeling formalism

Boolean networks consist of nodes with binary states that can be seen as on/off switches, and logical rules that determine state transitions. Each node is regulated by a defined set of upstream nodes and an associated update rule formalized through logical operators, such as AND, OR, NOT or XOR. The global state of the Boolean network is defined by the concurrent states of all nodes. Its transitions from one state to the next depend on the chosen semantics or update scheme, which determines what variables will have their state updated, and the functions they have been assigned. Two main semantics are typically used: synchronous updating, in which all variables are updated simultaneously for every transition, or asynchronous updating, in which subsets of nodes are updated at each step. In this work, we employed a synchronous update scheme to ensure deterministic trajectories and facilitate comparison with discretized time-series experimental data. Starting from an initial state, the Boolean network follows a trajectory that will eventually reach an asymptotic regime, or attractor, which can be a single steady state or a cycle pattern of states ⁶⁷.

Boolean modeling formalism is ideally suited for prediction tasks. Its inherent mechanistic information allows for hypotheses beyond the initial data scope. Additionally, its binary variables facilitate system perturbation and output interpretation. In practice, one can force the initial or constant state of one or more network nodes to either O/off or 1/on and observe how the other nodes are affected, i.e. which nodes remain stable and which ones are up- or downregulated.

Boolean naïve network construction

To build relevant Boolean networks of biological systems, a classical approach would rely on the manual curation of known interactions and their conversion into logical functions. To allow for the construction of larger networks and reduce our dependency on expert knowledge about the mechanisms described, we leveraged existing interaction information from curated signaling pathway databases. We used the Living Map app on the 3DEXPERIENCE platform (BIOVIA, Dassault Systèmes) to aggregate signaling pathways from the KEGG database ⁶⁸ and convert molecular interaction information into a single executable Boolean network. 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 ⁶⁹ thereby restoring connectivity while preserving directionality of regulation. To enhance immediate readability and interpretability, we connected the resulting network with the corresponding cellular population abundances analyzed by cytometry in the samples. We thus extended the idea of matching expression profiles of master regulatory genes with cell types from⁵⁰. These profiles being binary, they can be translated into Boolean functions, i.e. integrated directly into the network. We identified specific markers for each cellular population using the CellMarker⁷⁰ and Human Protein Atlas ⁷¹ databases and linked them to the rest of the network through a disjunctive function of their markers. The resulting networks could be referred to as naïve, as they have been assembled from aggregated sources of knowledge generated in various experimental contexts and species, with different measurement technologies.

Modeled experimental data

The naïve Boolean networks describing vaccine-induced immune responses were calibrated using longitudinal in vivo datasets generated in NHPs. For the MVA-induced immune responses, we used data from ⁴⁷. In this preclinical study, blood samples were collected from three cynomolgus macaques (Macaca fascicularis) at seven time points ranging from baseline (pre-vaccination) to one week after intradermal immunization. The sampling schedule included baseline, 3h, 6h, 24h, 48h, 72h, and 168h postadministration. Flow cytometry was performed to quantify immune cellular abundances and microarraybased transcriptomics to measure gene expression profiles on peripheral blood samples. On the other hand, the data used to calibrate the YF-17D vaccine-induced immune response model came from a separate cynomolgus macaque study described in ⁷². This dataset comprised blood samples from 12 animals collected at six time points: baseline and 1, 3, 7, 14, and 28 days after subcutaneous immunization. Samples were analyzed for immune cellular population abundances using mass cytometry, and gene expression levels were assessed using bulk RNA sequencing (RNA-Seq) transcriptomics. These two datasets provided time-resolved measurements of cellular and transcriptional responses, enabling the calibration of Boolean network dynamics to experimentally observed immune trajectories.

Context-specific calibration

To fine-tune naïve Boolean networks to the vaccine-specific contexts, we aligned them with experimental cellular abundance and gene expression data. We implemented a processing pipeline for these experimental data to suit the requirements of the Boolean formalism, converting continuous measurements into binarized time-series profiles. We filtered out the insufficiently strong signals, set a minimal fold change of 20% compared to baseline for at least one time point, and performed a Short Time-series Expression Miner (STEM) statistical filtering ⁷³ with a threshold of 0.05 on the gene expression data. As cellular abundance data contained neither enough animals for classical statistical testing nor sufficient dimensionality for the STEM algorithm, they were manually evaluated and subjected to a classical 2-fold change filter. All data were binarized using a 2-means clustering strategy ⁷⁴ on each expression/abundance profile. The resulting binarized data were then fed into the calibration algorithm named ZhegAlCal and presented in ⁷⁵ along with the naïve Boolean networks. The ZhegAlCal tool was parametrized to look for the network that maximized the overlap between the successive states of its nodes and the provided binarized profiles, while allowing for at most one intermediate network state between time points and imposing a fixed-point attractor.

Boolean network analysis, target identification and perturbation

Boolean networks can be analyzed both as through the perspective of graphs theory or dynamic simulation. We therefore combined topological and dynamic analyses to identify regulatory nodes with strong system-wide impact and to evaluate candidate intervention targets.

To quantify the size and connectivity of Boolean networks, we measured their number of nodes, interactions, and conjunctive clauses in the network regulatory functions. As a first attempt to identify their dynamically relevant nodes, we extracted their hubs, as defined in ⁷⁶. To refine this preliminary identification, we measured their Vertex Betweenness and Determinative Power, highlighted in ⁷⁷ as being associated with a strong dynamic impact. To bridge the gap with the actual networks’ dynamics, we computed the activity and effective graphs introduced in ⁷⁸ and ⁷⁹. Both concepts are presented as hybrid static/dynamic representations: the activity from one node to another is the probability that a perturbation of the former causes a perturbation in the latter in the next state. Quite similarly, the effectivity from one node to another captures the extent to which an incoming edge from the regulatory node is, on average, necessary to determine the value of the regulated node. Interestingly, this second notion has been recognized as an alternative definition for collective canalization, itself a generalization of the concept of canalization introduced in ⁸⁰ and shown to be inherent to a vast majority of biological processes ⁸¹. Based on these graphs, we implemented a search and display of the nodes that combined large numbers of output edges and high average activity/effectiveness, which we posited as being the variables with the strongest dynamic impact on their downstream nodes. We pushed the concept further as we then computed an algorithm to iteratively find the nodes with a large activity/effectiveness upstream of the identified nodes to define highly active or effective signal paths. These flows could indicate preferential targets to strongly affect the dynamics of a Boolean network.

Finally, we designed a workflow to perturb the Boolean networks by forcing a constant activation (on state) or inhibition (off state) onto a selected set of nodes suggested by the static analyses and comparison, and visualize the resulting changes in their dynamics. We additionally favored actionable target nodes through MVA genetic alterations. The changes are displayed as a downregulated vs. stable vs. upregulated heatmap from unperturbed to perturbed conditions. However, such a heatmap can be investigated in detail only for a limited number of nodes. To extend our analysis to a system-wide scale, we adapted Gene Ontology (GO) enrichment analysis ⁸² to our work. GO enrichment analysis is generally used to identify overrepresented terms in a gene set compared to a reference set, linking them to biological processes, molecular functions, and cellular components. For a Boolean network in a given state, we can extract the upregulated genes that went from off at baseline to on, or the downregulated genes that went from on to off and generate the resulting enriched GO terms; compared to all nodes in the network. In this work, we concentrated our analysis on biological processes using the 2023 Gene Ontology database ⁸³ via the GSEApy Python library ⁸⁴.

Results

The construction of Boolean models of vaccine-induced immune responses comprised several sequential steps. First, naïve Boolean networks were generated by converting and merging several signaling pathways into Boolean logic and enriching the resulting network with oriented protein-protein interactions and immune cellular populations. Context-specific experimental data were then provided as binarized time series and used to calibrate the naïve Boolean networks. The calibrated Boolean networks were subsequently simulated and analyzed to produce mechanistic insights and to identify and evaluate candidate perturbation targets. The complete workflow is described in Figure 1.

Figure 1:

Construction of a naïve Boolean network of the immune response to MVA

Although its cellular receptor remains unknown, MVA enters the cells similarly as VACV by direct fusion with the plasma membrane through a fusion complex or macropinocytosis ^85,86, and preferentially infects antigen-presenting cells (APCs), such as DCs and macrophages ^87–89. Upon infection, DCs can follow two main fates: activation through the cGAS/STING-mediated cytosolic DNA-sensing pathway, as shown in a murine model in ⁹⁰ and confirmed in human DCs in ⁹¹, or cell death via apoptotic mechanisms. In addition, non-infected DCs can be indirectly activated through cytokine production by neighboring infected DCs ⁹¹. In murine and human macrophages, MVA infection triggers diverse apoptotic pathways ^92,93 and activates the NF-κB pathway ^92,94 Interestingly, the parental VACV strains Copenhagen ⁹⁵ and Western Reserve (WR) ⁹⁶ encode several intracellular immunomodulatory proteins that inhibit NF-κB signaling, such as A52 ⁹⁷ and N1 ⁹⁸, which are not conserved in the MVA genome. This difference could contribute to the observation that MVA, but not WR, activates NF-κB in human macrophages and fibroblasts ⁹². This previous knowledge was used to define the backbone of a representative naïve Boolean network, i.e. its topology or which nodes compose the network and how they interact. We merged the Cytosolic DNA-sensing (which includes type I interferon (IFN) signaling), apoptosis, and NF-κB signaling pathways from KEGG ⁶⁸ into a single gene regulatory network using BIOVIA Living Map, a software on the Dassault Systèmes 3DEXPERIENCE platform (BIOVIA, Dassault Systèmes). These pathways encompass key mechanisms important in how MVA infects and is sensed by host cells, and thereby induces essential innate immune responses. In addition, NF-κB and type I IFN signaling also play a crucial role in amplifying these inflammatory and antiviral responses in both infected and uninfected (bystander) cells.

Network calibration with experimental datafrom a NHP immunogenicity study

To calibrate the naïve network, we used data from the preclinical study described in ⁴⁷. Healthy adult male cynomolgus macaques (Macaca fascicularis) were immunized intradermally with a recombinant MVA expressing eGFP. The experimental data underwent preprocessing. In particular, we implemented a data filtering strategy compatible with datasets derived from small numbers of individuals. The resulting binarized time-series data were used as input to the ZhegAlCal calibration algorithm ⁷⁵ together with the naïve Boolean network (Figure 1). ZhegAlCal was to identify the network configuration that maximized the overlap between successive node states and the corresponding binarized experimental profiles. Calibration allowed for at most one intermediate network state between consecutive experimental time points and imposed convergence toward a fixed-point attractor. The calibration Python script was run in less than one day on a Windows 10 workstation equipped with an 11th Gen Intel Core i7-11850H CPU (2.50 GHz) and 32 GB RAM.

Characterization of the calibrated Boolean network of the immune response to MVA

The calibrated Boolean network describing the innate immune response to the MVA vaccine in NHPs comprised 211 nodes, including 200 genes, proteins, or molecules from the Cytosolic DNA-sensing, apoptosis, and NF-κB signaling KEGG pathways, as well as 11 immune cellular populations. Nodes were regulated by up to 11 inputs. The network is illustrated in Supplementary Figure S1 and is available for simulation, perturbation, and analysis on the CellCollective platform. Each cell type was connected to the gene regulatory components of the network through logical functions based on gene markers obtained from CellMarker⁷⁰ and the Human Protein Atlas ⁷¹.

To ensure that the calibrated Boolean network covered a broad spectrum of immune processes, we performed a static Gene Ontology (GO) enrichment analysis of the 200 genes derived from the three merged selected KEGG pathways using the human genome as reference (due to the lack of a cynomolgus macaque GO reference) (Figure 2). We observed a comprehensive overview of biological processes classically involved in responses to viral infection⁹⁹. As expected, these included processes intrinsic to the NF-κB and apoptosis pathways. Importantly, enrichment also encompassed cytokine-related processes, particularly TNF signaling and production, as well as a broader inflammatory processes and alternative sensing pathways, supporting the ability of the calibrated network to capture systemic immune responses beyond infected APCs.

Figure 2:

The dynamics of the simulated Boolean network for immune cellular populations are shown in Figure 3a. We employed synchronous update semantics to ensure deterministic trajectories between any successive states, which also tend to yield shorter and more interpretable trajectories while preserving essential information about the network ¹⁰⁰. During calibration, intermediate states were introduced between certain experimental time points (between 24h-48h, 48h-72h, and 72h-168h) to maximize the overlap with the provided data. To evaluate the efficiency of the calibration process, faithfulness with the calibrated data was assessed by comparing the simulated profiles of the network nodes that appear in this data with their binarized experimental values (Figure 3b). We obtained satisfactory average ratios of correctly predicted states of 87% for the differentially expressed genes (not shown) and 86% for the differentially abundant cells (Figure 3c). To further evaluate the biological relevance and the comprehensiveness of the simulated trajectory of the Boolean network, we computed a GO enrichment analysis for the upregulated nodes for each one of its successive states after baseline. The enrichment results on the successive time points showed a broad response to the vaccine administration (Supplementary Figure S2). The GO terms that were upregulated after one update of the network covered a variety of processes without any clear trend. However, the next state, corresponding to 6h post-administration, highlighted a substantial activation of cellular migration processes, including chemotaxis, consistent with immune cell recruitment from the bone marrow to the vaccine injection site through the bloodstream, which is part of the inflammatory response ¹⁰¹. The subsequent state, associated with the 24h time point, showed strong enrichment of T cell-related processes, particularly T helper 1 responses, along with regulation of associated cytokines, such as IL-2, in agreement with observations reported in ⁴⁷. After two additional updates, the Boolean network state reflected regulation of both innate and adaptive immune mechanisms, with prominent upregulation of host responses against viral and double-stranded RNA (dsRNA) stimuli. This aligns with dsRNA production during later phases of poxvirus replication ¹⁰², and with reports that a modified MVA strain producing excess dsRNA early during infection induces stronger innate immune response ¹⁰³. The last four states of the trajectory maintained strong antiviral signatures and showed enrichment of hydrogen peroxide-related processes up to one week post-infection, suggesting sustained release of reactive oxygen species (ROS) to fight the infection. In contrast, decrease ROS production has been reported as an immune evasion mechanism of VACV ^104,105.

Figure 3:

Overall, analysis of the unperturbed network dynamics confirmed that the calibrated Boolean network recapitulates the experimentally observed sequence of immune events induced by a parental MVA in the reference study ⁴⁷ and in related contexts ¹⁰⁶. These results support the ability of the Boolean formalism to qualitatively capture key features of the complex early vaccine-induced immune responses, in general and to MVA specifically as proof-of-concept, despite its inherent simplifications.

Challenging the Boolean network with vaccine mutations

To further validate the Boolean network of immune response to MVA and support its biological relevance and legitimize its use, we evaluated its predictive capacity under previously unseen conditions. We perturbed the system to mimic immunization with modified MVA vectors described in the literature to be more immunogenic than the parental MVA used to construct the model. The outcomes of these in silica perturbations were then compared with the corresponding experimental data.

The Boolean network we built and calibrated for our application is a model of the host immune responses and exclusively integrates host genes, proteins, and cellular populations, corresponding to an averaged representation of three animals. To mimic the alterations of the MVA genome that constitute our experimental references for validating the heterologous predictions, we needed to identify their major interactors within the host proteome in existing literature. Several mutant MVA vaccines with deletions in immunomodulatory genes, including A40R, A46R, C6L, K7R, and N2L, have been described and evaluated for replication and immunogenicity ^{23,25,26,28,29}. The host targets of some of these viral proteins are well characterized. The C6 protein interacts with the TBK-1 adaptor proteins TANK, TBKBP1, and AZI2, thereby inhibiting IRF3 and IRF7 activation at the level of MAVS, IKBKE, and TBK1 ¹⁰⁷. The K7 protein also inhibits IRF3 and IRF7 activation via the IKBKE/TBK1 axis by binding DDX3X host protein and also suppresses TLR-induced NF-κB activation by its targeting of IRAK2 and TRAF6 ¹⁰⁸. Finally, the A46 protein interferes with IRFs, NF-κB, and MAPKs activation by sequestering adaptor proteins MYD88, TIRAP, TICAM1, and TICAM2 ¹⁰⁹, thereby preventing the assembly of MYD88 and TIRAP and stopping their signaling cascade ¹¹⁰. All of the aforementioned MVA-host interactions are inhibitory. Thus, deleting the given viral gene is expected to reduce the inhibitory effect of the corresponding viral protein. We mimicked this process by forcing the activation of target host genes in the Boolean network.

We next evaluated in silica immune responses to the triple mutant MVA AC6L/AK7R/AA46R and compared them with experimental data from in vitro and in vivo mouse studies using two recombinant MVA AC6L/AK7R/AA46R expressing Chikungunya and Ebola virus antigens ^23,28. As readouts, we selected genes whose expression we had experimentally quantified in comparison to the unmutated parental MVA, including IFN-β (IFNB1), IL-6, TNF-γ (TNF), and CCL5 (RANTES) inflammatory cytokines and chemokines, as well as two IFN-stimulated genes (ISGs): the type I IFN-inducing cytosolic pattern recognition receptor RIG-I, and the chemokine IFN-γ-induced protein 10 IP-10 (CXCL10). We displayed the changes between the wild-type and triple-mutant MVA perturbed model dynamics in a downregulated (i.e., that went from 1 to 0) vs. stable vs. upregulated (that went from 0 to 1) heatmap for these selected genes and the network cellular populations (Figure 4a). The binary trajectories underlying this confrontation are given in Supplementary Figure S3. The perturbed model was obtained by forcing the activation of the following genes: TICAM1, TICAM2, TBK1, IKBKE, MAVS, TRAF6, and the MYD88:TIRAP complex, corresponding to the deletion of the three viral genes, all known to show inhibitory effects.

Figure 4:

We observed that the triple mutant MVA ΔC6L/ΔK7R/AA46R upregulated the expressions of IFN-β and of the proinflammatory cytokines IL-6 and TNF-α, as well as IP-10, consistent with enhanced inflammatory and innate responses. Cellular populations showed less drastic perturbations of their abundance. We saw an upregulation at earlier time points for plasmacytoid DCs (pDCs) and B cells, and at later time points for natural killer (NK) cells. Towards the end of the simulated trajectories, B cells and classical monocytes showed a downregulation, whereas MDSCs and granulocytes showed less intelligible fluctuations. Overall, these observations were concordant with the reported experimental phenotype of the mutant MVA-induced responses assessed in vitro using the human monocytic cell line THP-1, and in vivo in mice ^23,28, particularly the early upregulation of IFN-β, TNF-α, and IL-6. The lack of perturbation we observed for RIG-I also aligned with experimental observations, unlike some of the changes we saw for IP-10 and, conversely, did not see for CCL5 (RANTES). It suggested that the Boolean network gave a credible description of the molecular and cellular innate mechanisms involved in these responses. To widen our investigation and look at the full Boolean network dynamics, we also displayed the system-wide downregulated vs. stable vs. upregulated heatmap, with all network nodes making up the y-axis (Figure 4b). This illustration highlighted the fact that few perturbations (in this case just 7) strongly affected the network behavior, in both permanent and transitory ways.

To provide a biological interpretation for these perturbations, we computed the GO enrichment analyses for the upregulated genes at each time point of the mutant MVA-induced dynamics (Supplementary Figure S4). By comparing, update-by-update, these GO enrichment results with the enrichments of the parental MVA-induced Boolean trajectory (Supplementary Figure S2), we noted a similar cascade of early events followed by a resembling, yet presumably broader (less T helper-oriented), T cell response, a longer-sustained upregulated acute inflammation, and an earlier antiviral response at 24h postadministration. The network states that led to the 48h time point translated a continuation of these features, thereby converging back to the unperturbed response. We could, however, observe an earlier upregulation of ROS secretion. Finally, the later network states in perturbed conditions closely resembled the unaltered dynamics despite a notably stronger downregulation of apoptotic processes.

Together, these network perturbation analyses demonstrated the predictive capacity of the Boolean model for MVA-induced immune responses, and support the validity of the modeling strategy. Next, we used the Boolean model of immune response to MVA to evaluate novel MVA mutants derived from expert knowledge.

Guiding the design of new mutant MVAs using the Boolean network

We next harnessed the Boolean model of the MVA-induced response to virtually evaluate two new mutant MVA vectors. VACV N2L (MVA021L) has been reported to encode a protein that inhibits IRF3 activation downstream of TBK1 after IRF3 translocation into the nucleus ¹¹¹, and deletion of N2L from MVA has been associated with increased immunogenicity in mice ²⁶. Using the Boolean model of MVA-induced response, we tested the effects of combining N2L deletion with deletions of C6L, K7R, and A46R by forcing IRF3 activation, in addition to activation of TICAM1, TICAM2, TBK1, IKBKE, MAVS, TRAF6, and the MYD88:TIRAP complex. Consistent with in vitro data obtained in THP-1 cells and monocyte-derived DCs ²⁶, we observed a minor and transient change in IFNB1 behavior for the single N2L mutant compared to the wild-type MVA, but not of TNF. However, there was no effect on immune cell subset abundances (Figure 5a), nor was there a difference between the quadruple MVA ΔN2L/ΔC6L/ΔK7R/ΔA46R mutant and the triple MVA ΔC6L/ΔK7R/ΔA46R mutants when compared to the parental MVA (Figure 5b and Figure 4a).

Figure 5:

The MPXV OPG147 gene encodes a protein that has been shown to inhibit the cGAS-MITA/STING-mediated innate immunity by preventing its ISGylation, oligomerization and trafficking ¹¹². This gene is conserved among poxviruses; VACWR140 (A21L) and MVA 131L are orthologs ²¹. Substitution of three alanines at positions 55, 116, and 117 of VACV A21L gene has been reported to maintain VACV competence for viral entry and virus-induced cell-to-cell fusion, while abrogating its interaction with MITA/STING1, resulting in enhanced innate immune responses compared with wild-type VACV in vitro (THP-1 cells) and in vivo in mice ¹¹². To model this mutation, we activated MITA/STING1 in the Boolean model of the response to MVA to mimic immunization with an MVA A21L-F55A/T116A/T117A mutant. This perturbation alone led to transient IFNB 1 upregulation (like in mice) and increased pDC abundance compared with parental MVA (Figure 5c). However, the quadruple MVA A21L-F55A/T116A/T117A/ΔC6L/ΔK7R/ΔA46R mutant and the triple MVA ΔC6L/ΔK7R/ΔA46R mutant when compared to parental MVA did not differ qualitatively (Figure 5d and Figure 4a).

No qualitative differences were observed in silico between the triple MVA AC6L/AK7R/AA46R mutant and the two new quadruple MVA mutants (MVA ΔN2L/ΔC6L/ΔK7R/ΔA46R and MVA A21L-F55A/T116A/T117A/ΔC6L/ΔK7R/ΔA46R). However, the inherent limitations of Boolean modeling, which captures only binary state changes, prevent the capture of quantitative differences between them that could result from additive or synergistic effects of N2L deletion or alanine substitutions in A21L in combination with deletions of C6L, K7R and A46R. Given that not only the quality and dynamics of innate immune responses, but also their magnitude, modulate adaptive immune responses, the construction and in vitro and in vivo testing of these two new MVA quadruple mutants, identified using the Boolean model of response to MVA in a supervised, expert-guided manner, might still be worthwhile. Next, we investigated the potential of Boolean modeling for unsupervised rational discovery of additional MVA variants.

Improving MVA immunogenicity through comparison with the YF-17D-induced immune response

Construction of a consensus naïve network

To fully exploit the model of responses to MVA as a discovery tool to improve its immunogenicity, we compared it with a model of responses to YF-17D, one of the most successful vaccines developed to date. This comparison first required the construction of a shared naïve network between MVA and YF-17D.

The YF-17D vaccine has been shown to infect several DC subtypes, such as myeloid DCs (mDCs) and pDCs ^38–44 Upon infection, these DCs signal through multiple Toll-like and RIG-I-like receptors⁴¹. Immature DCs undergo apoptosis upon infection, unlike mature DCs ³⁸. Notably, the production of IFN-a by pDCs during YF-17D infection requires activation of the mTOR pathway, which is also necessary for optimal antigen-specific CD8 T-cell responses ¹¹³. Based on this literature-derived knowledge, we enriched the initial MVA-specific naïve Boolean network (constructed from 3 KEGG pathways) with the Toll-like receptor, RIG-I-like receptor, and mTOR signaling pathways, also obtained from KEGG ⁶⁸ and imported into BIOVIA Living Map (BIOVIA, Dassault Systèmes). The resulting consensus naïve network recapitulated the main biological entities and mechanisms that characterize both MVA and YF-17D vaccine-induced innate immune responses.

Separate calibrations

Using this shared topology, MVA and YF17D vaccines-specific immune responses were defined through separate calibrations with dedicated preclinical datasets. For MVA-induced responses we used the gene expressions and cellular abundance data from ⁴⁷, whereas for YF-17D-induced responses we used the gene expressions and cellular abundance data from a separate cynomolgus macaque study from ⁷². To maximize comparability, both networks were calibrated using only data from overlapping time points: baseline, 1 day, 3 days, and 1-week post-injection. The ZhegAlCal tool from⁷⁵ was applied to the consensus naïve network and separately with the processed and binarized MVA and YF-17D datasets. It was parametrized to allow for at most one intermediate network state between experimental time points and to assume convergence to a steady-state attractor, consistent with the studied time scale. The respective Python calibration scripts were run in 12 and 37 hours for the MVA and YF-17D networks, respectively, on a Windows 10 workstation equipped with an 11th Gen Intel Core Î7-11850H CPU (2.50 GHz) and 32 GB RAM. Both calibrated Boolean networks are provided in Supplementary Figure S5 and are available on the CellCollective platform.

Comparative dynamics analysis

To compare the behavior of the MVA and YF-17D vaccine-induced response Boolean networks, we computed their dynamics and displayed binary trajectories of cellular abundances (Figure 6a and Figure 6d), the corresponding experimental binary profiles (Figures 6b and 6e), and their comparison (Figure 6c and Figure 6f). We saw that the calibration algorithm exploited the potential additional network update between experimental time points to better match the provided data. The MVA network included one additional state between 3- and 7-days post-injection, whereas the YF-17D network included two intermediate states, between 1 and 3 days and between 3 days and 1-week post-injection. The two models accurately recapitulated the experimental data used for calibration, correctly predicting the states of 92% and 86% of differentially expressed genes in the MVA and YF-17D networks, respectively (not shown), as well as 98% and 100% of differentially abundant cellular populations (Figure 6c and Figure 6f).

Figure 6:

We performed a GO enrichment analysis for both vaccine-induced response networks at every state after baseline (Supplementary Figure S6). The successive states of the MVA response network translated an early (1-3 days post-injection) upregulation of calcium-related processes, STAT phosphorylation, and global downregulation of leucocyte apoptosis. The GO terms suggested a concentrated cytokine and chemokine environment that induced cellular migration, along with early regulation of T lymphocyte activation and antiviral processes. The later states still showed sustained T cell activation, cytokine and chemokine-related as well as antiviral processes. Yet, the response seemed to have become broader, with a transition toward adaptive immune mechanisms, including translated B cell differentiation, activation, and NK cell activity. The final state (1-week post-injection) expressed this broad response with a weaker antiviral response but enhanced IL-2 production. Overall, when enriched with three additional signaling pathways to cover a broader spectrum of MVA-inducible mechanisms (but still calibrated with the data from ⁴⁷), the unperturbed response retained its properties. The GO terms associated with the trajectory followed by the YF-17D response network differed from that of MVA, as expected from experimental data ³⁷. They showed strong metabolic activity and signaling through the NF-κB signaling pathway at 1 day post-injection. The following state revealed a broader response where the upregulated metabolic processes persisted alongside an enhanced T cell activity and PI3K-Akt signaling, as well as downregulated apoptotic processes. We also noticed the contribution of a type II IFN-mediated response. At 3 days post-injection, apoptotic processes remained downregulated. In addition, the GO terms indicated an important role for Fc receptor-induced signaling. The last two states translated substantial signaling through the MAPK pathway, intertwined with an enhanced response to stress. A variety of other phenomena were also observed one week after immunization with the YF-17D vaccine, including cellular antiviral processes, granulocyte chemotaxis, and a response to oxygen species that suggests ROS secretion.

Comparative static analysis

We next performed a static analysis of the two Boolean networks. Because the MVA and YF-17D datasets included different sets of filtered cellular populations, the MVA-induced response Boolean network had one additional cell-type node compared with the YF-17D network (Supplementary Figure S5) (329 vs. 328 total nodes). The MVA network also had 2 more interactions and 5 more conjunctive clauses that connected these cellular populations to the rest of the nodes. Interestingly, this difference did not extend to the full calibrated networks. Indeed, the YF-17D-induced response network exhibited a larger number of interactions and conjunctive clauses (583 and 509, respectively) compared with the MVA network (576 and 481, respectively).

Moving on to the graph theory measures, both networks had the same hub nodes (as defined in⁷⁶): FOS, GZMB, IRF3, IRF5, MTOR, NFKB1, NFKB1:RELA, TP53, TRAF2, and TRAF6. The refined Determinative Power (DP) and Vertex Betweenness (VB) metrics were confronted for the MVA- and YF-17D-induced responses network (Figure 7a and Figure 7b). We annotated the dots corresponding to the nodes that displayed high VB and DP. Both networks also attributed high VB and DP to the same node, AKT3. In addition, this analysis highlighted the shared importance of the RELB:NFKB2 complex and the INSR, TNFRSF1A, RPS6KB1, PIK3CA, and MAP3K7 nodes. We computed the effective graphs ⁷⁹ of the MVA and YF-17D-induced response Boolean networks. We then displayed the average effectiveness of their nodes as a function of the number of output edges (Figure 7c and Figure 7d). We annotated the dots that combined large numbers of output edges and high average effectiveness. Once again, the highlighted nodes for both MVA and YF17D vaccine-induced networks aligned with their hubs, with the additional node RELB:NFKB2. We also noted the presence of CD14 among them for the MVA-induced response, which could be a bias due to the large number of monocyte populations that are connected to this marker. Finally, as in ⁷⁹, we measured the average effectiveness of the regulatory functions in both Boolean networks as a function of their number of input variables (Figure 7e). We obtained values that aligned with what is presented in ⁷⁹ for Boolean networks stored on the CellCollective repository, and that were substantially smaller than for randomly generated networks. According to the authors, this reduced effectivity implies collective canalization, which is an indicator of robustness to perturbations and a property shared by most real biological systems.

Figure 7:

MVA-induced immune response improvement via Boolean network comparison

Both the cellular population dynamics (Figure 6a and Figure 6b) and the GO enrichment analyses of system-wide network states (Supplementary Figure S6) showed major differences in how the hosts (cynomolgus macaques here) respond to the MVA and YF-17D vaccines, as expected. An overview of these differences indicated a less extensive, or potentially delayed, adaptive response in favor of innate immune mechanisms and stress-induced MAPK signaling for YF-17D, as previously described ³⁷. Yet, vaccine optimization has proven to be a subtle procedure, and the dissimilarity between the modeled responses was not specific and precise enough to point out a clear improvement strategy straightforwardly. Next, we investigated which components of this complex system could best be affected to orient the response in the desired direction, and how to do so.

Building on the identification of dynamically impactful network nodes, we developed an automated tool to find highly effective signal paths leading to these nodes. Nodes within these paths were conjectured to influence the network dynamics when targeted, as their perturbation significantly affects the state of their downstream nodes toward one or more of the identified dynamically impactful nodes, which, in turn, powerfully impact the entire network. We computed the highly effective paths leading to the nodes of interest identified and annotated in Figure 7c (Figure 7f). Historically, two strategies have been used to improve MVA viral vaccine vectors: the insertion of exogenous co-stimulatory molecules into the MVA genome and the deletion of some of its remaining immunomodulatory genes²². The overlap between tested co-stimulatory molecules, cytokines, and chemokines (IL-2, IL-12, and IL-15, IFN-γ, OX40/OX40L, a triad of B7-1, ICAM-1, and LFA-3 molecules, a triad of CD80, CD86, and CD83 molecules, CD40L, or GM-CSF) and the nodes included in at least one effective path yielded no results. On the other hand, many of the tested poxvirus deletions have been shown to affect host genes within these paths. After filtering out the deletions concerning viral genes that are not conserved in the MVA genome (which differs distinctly from the VACV it originated from ²¹), we intersected the host interactors of the remaining viral genes with the nodes in the highly effective paths. The matching nodes in the host response network are CHUK, IKBKB, MAVS, NFKBIA, TRAF6, and ZBP1. These are antagonized by proteins encoded by B14R, C6L, K1L, K7R, and E3L viral genes, respectively. In addition, the viral poxvirus gene F17R, which is conserved in the MVA genome ²¹, encodes the F17 protein recently shown to impair MTOR signaling upon infection ¹¹⁴. On the other hand, despite being (at least partially) conserved in the MVA genome, the functional integrity of the K1L gene, and therefore the effect of the K1 protein on the NFKBIA host gene, is unknown. We thus simulated the immunization with mutant MVA viral vectors with deletions in B14R, C6L, K7R, E3L, or F17R, which have never been considered, designed, constructed and evaluated on their own so far, to our knowledge. As all five viral proteins are inhibitors of the associated processes, we mimicked their deletion by a forced activation of their interactors. To evaluate the changes in the resulting network dynamics, we computed their downregulated vs. stable vs. upregulated heatmap between unperturbed (parental MVA) and perturbed (new mutants) conditions. In addition to the cellular populations, we highlighted the perturbations of genes involved in the induction or function of the stress-related JNK and p38 MAPK pathways, which are involved in the response to YF-17D. These include the IL-α (IL1B), TNF-α (TNF), and their receptors (IL1R1, TNFRSF1A), as well as the MAPK kinases MKK3 (MAP2K3) and MKK7 (MAP2K7), the MAPK kinases kinases MEKK1 (MAP3K1), ASK1 (MAP3K5), and TAK1 (MAP3K7), and finally the actual JNK1 (MAPK8) and p38α (MAPK14) genes. We also added the SIN1 (MAPKAP1) gene that is part of the mTORC2 complex, thereby bridging the gap between the stress-related MAPK and mTOR pathways, with the latter being an essential component of the YF-17D vaccine-induced response ¹¹³. We observed no noticeable changes for the ΔB14R, ΔC6L, and ΔE3L simulations (data not shown), in contrast to ΔK7R and ΔF17R MVA mutants’ trajectories (Figure 8a and Figure 8b). The binary trajectories underlying these latter comparisons are given in Supplementary Figure S7. The perturbations corresponding to the K7R deletion in MVA induced a slight decrease in NK and CD8 T cells at 1week postinjection. Interestingly, we also observed an increase in IL-1 β (logically followed by its receptor) and MAP3K1 that could respectively initiate and contribute to the JNK MAPK pathway (KEGG PATHWAY: MAPK signaling pathway - Homo sapiens). On the other hand, the perturbations corresponding to the F17R deletion in MVA provoked an upregulation of MAPKAP1 at later time points. In conclusion, we have demonstrated that our model construction and analysis strategy is a powerful tool for identifying new vaccine candidates, here two new recombinant MVAs that warrant further characterization, on their own as single mutant or combined with other deleted viral genes.

Figure 8:

Discussion

We previously designed and implemented a methodology to build and analyze large Boolean networks ⁷⁵. We merge signaling pathways to form a naïve network with up to several hundreds of nodes and use a calibration tool that successfully aligns such a naïve network with context-specific experimental data. Here, we applied this methodology to provide mechanistic insights and vaccine improvement strategies with the virtual design and evaluation of new vaccine candidates. We used a preclinical trial with MVA as proof of concept. We first validated the methodology against experimental observations. The calibrated MVA-induced response network recapitulated experimental characteristics. We then validated the model and the modeling methodology through network perturbation. Its predictive capabilities were confirmed, as its simulation recapitulated the immunogenicity of a previously-reported mutant MVA vaccine, without the corresponding experimental data being used to construct the model. We finally exploited networks to uncover new vaccine strategies. New candidate mutant MVA vaccines were discovered by a rational design based on Boolean models and evaluated virtually. We tested deletions or point mutations of viral immunomodulatory genes derived from two sources: previous knowledge of host/poxviruses interactions on the one hand, and mechanistic insights gained from the analysis of the models on the other. We successfully perturbed the MVA input signal to mimic the YF-17D-induced immune response by identifying MVA-induced signaling paths with high dynamic impact. Therefore, Boolean modeling of immune responses to vaccines appears to be a useful and powerful tool for accelerating vaccine development by enabling the discovery, tailoring, testing, and ranking of new vaccine candidates to be further evaluated in vitro and in vivo.

Keeping in mind that the choice of the pathways that compose the naïve network is of paramount importance, as the biological entities they contain should allow for an exhaustive yet adequately specific description of the underlying phenomena, we built our initial Boolean model of the innate immune response to the MVA vaccine around three signaling pathways, mainly involved in detecting and initiating the immune response to viral stimuli. Nevertheless, the scope of the integrated KEGG NF-κB pathway simultaneously provides the necessary elements to depict the various immune, particularly inflammatory, processes occurring beyond infected cells. Overall, the static GO enrichment performed on its nodes showed a satisfactory balance between specificity and exhaustivity for describing MVA-induced immune responses. Moreover, the calibration process successfully overlapped its trajectory with longitudinal binarized experimental measurements. The comprehensive structure combined with the efficient calibration yielded a Boolean network whose successive states represented a series of biologically relevant and experimentally assessed processes¹⁰⁶.

Simulating this network in perturbed conditions induced systemic modifications that strongly affected the network, and specifically the selected output markers (immune cell subsets abundance). The triple mutant MVA ΔC6L/ΔK7R/ΔA46R, modelled through the forced activation of the few host genes whose proteins have been reported to be inhibited by the C6, K7 and A46 viral proteins, notably induced changes over the entire network. These perturbations highlight an overall broader T-cell response, an upregulated acute inflammation, earlier upregulation of virus-induced innate immune responses and ROS secretion, and a downregulation of apoptosis-related processes at later time points, compared to the parental MVA. Such enhanced innate immune response that manifested through these virtual perturbations, particularly the early upregulation of IFN-β, TNF-α, and IL6, matched experimental data from ^23,28.

Finally, we added three more KEGG pathways to the naïve network to obtain a consensus basis for both the MVA and YF-17D vaccines-induced responses Boolean networks. Together, the six pathways provided the necessary information for describing the variety of processes involved in the immune response to these live-attenuated vaccines. The static GO enrichment on the gene regulatory network resulting from the merger of the six pathways showed that these processes include mechanisms ensuing directly from the selected pathways, as expected, but also broader transversal phenomena such as inflammation and the production and response to a variety of cytokines, including TNF-α. An extensive analysis of the topological properties of both networks revealed a close resemblance. A slight difference, however, that could provide some insights into what causes the observed responses to differ, is the number of interactions, which is higher in the YF-17D-induced response network. This translated a stronger connectivity within the network, but could also imply a higher redundancy in the paths from one variable to its downstream nodes. Redundancy is generally associated with robustness but also with collective canalization ⁷⁹. Collective canalization implies that a subset of the input variables for every regulatory function governs its output. In⁷⁹, the authors also present the notion of effectiveness and associate it with collective canalization. We computed the effective graph for the MVA and YF-17D networks. We observed values that were, on average, substantially smaller than what is displayed by randomly designed networks, suggesting biologically relevant regulatory rules in both networks, as well as an overall tendency to robustly withstand perturbation. In this context, the numerous system-wide changes observed when forcing the activation of just a few nodes in the initial MVA Boolean network suggest that this model (and our underlying construction methodology) grasps the biological relevance of these immunomodulatory genes. Coming back to the MVA and YF-17D network comparison, we additionally note that the average effectivity was always smaller for the YF-17D-induced response network (except for 7-input functions), corroborating the slightly higher redundancy observed for this network.

The respective network trajectories, on the other hand, could hardly be more different. We observed almost opposite behaviors for several essential cellular populations present in both networks, such as (classical) monocytes, B cells, NK cells, CD4 and CD8 T cells. The successive GO enrichment analyses confirmed these observations, as the computed terms for all time points after baseline translated different biological processes. The successive MVA network states described a standard immune response, with early innate and inflammatory mechanisms that induced and guided an adaptive response before slowly fading away to leave only the adaptive processes. The trajectory followed by the YF-17D network, on the other hand, yielded less intuitive GO terms. Overall, we observed a response that was initiated by signaling through the NF-κB and Pl3K-Akt pathways, which then involved a variety of mostly innate immune mechanisms, even towards later time points and up to the final time point in this study at 1-week postinjection. Adaptive immunity processes remained notably scarce over the entire trajectory within this short timeframe. Interestingly, this observation overlaps with what is described in ³⁷ where the transcriptional signature of the yellow fever vaccine YF-17D was compared with the signatures of 12 other vaccines, including a poxvirus-based smallpox vaccine (VACV, not MVA) and a recombinant MVA against tuberculosis. They also exhibited a delayed immune response, with antiviral processesjust starting at 3 days post-injection and reaching their peak after 1 week, which was also the only time point where the differential enrichment of the inflammatory processes was statistically significant.

Despite the differences between the MVA and YF-17D-induced responses and corresponding networks, no clear optimization strategies for the MVA viral vector came to mind from the comparison. We looked for components we could efficiently affect. We computed the highly effective signal paths that gave us the nodes most likely to strongly impact the network dynamics. We overlapped these nodes with the host genes known to interact with MVA viral genes. Among these viral genes, we found B14R, C6L, K7R, E3L, and F17R to possess interactors within our highly effective paths. The first four are known and their single deletions from poxviral genomes have been tested with optimization of their immunogenicity in mind ¹¹⁵. The last one, F17R, has only relatively recently been linked with an inhibitory effect on the mTOR signaling pathway, which immediately stimulated our interest, given the preponderant role of this pathway in the YF-17D response. We thus computed the perturbed responses mimicking a deletion in each one of the genes. Whereas the ΔB14R, ΔC6L, and ΔE3L simulations yielded no sensible observations, the deletion of K7R, and F17R to a lesser extent, showed interesting consequences. We noted higher levels of IL-1β and MAP3K1, which is involved in the stress-induced JNK MAPK pathway, for the former and slightly higher MAPKAP1 for the latter, which suggests an enhanced cross-activity with the mTOR pathway.

In truth, several aspects of our methodology could be improved and should perfect its outcomes. An important limitation of this study might be the Boolean approximation itself. Evidently, not all biological processes, especially gene expressions and cellular abundances, can realistically be described as on/off switches. This strong approximation can dilute subtle changes in larger trends, which might be what prevented us from observing more differences between the different perturbation scenarios. More generally, differences in the magnitude of observed phenomena can hardly be distinguished. To still capture as much biologically relevant information as possible, the data processing and binarization steps are of paramount importance. Notably, the 2-means clustering we employed here implies that an identical expression value in two different conditions can be assigned to either cluster depending on the other values. Again, this could explain why such little changes occurred between some scenarios and why some observations slightly diverged from experimental reality. Substantial improvement could be provided by a shift from Boolean to ternary ¹¹⁶ or multivalued ^117,118 logic. Another refinement could be provided by changing not the number of potential states but the way they are updated. Examples include the addition of probabilistic activation and inhibition rates ¹¹⁹, or the application of the Most Permissive Boolean Network semantics ¹²⁰.

Another vital parameter in our approach is what we use to build the networks. First, the knowledge we use as the backbone of these networks is of paramount importance. Despite its satisfactory coverage of the biological phenomena involved in the innate immune response to the MVA vaccine, the choice of signaling pathways to compose the naïve network admittedly skews the final model and the conclusions that can be drawn from it. Indeed, the Boolean networks of the MVA-induced immune response, calibrated with the same dataset but based either on a 3 - or 6-pathway naïve network backbone yield similar yet not identical simulated behaviors. In addition, as we remain limited in the number of pathways to add to our network backbone (mainly for computational reasons), some processes, but also relevant markers such as the macrophage inflammatory proteins (MIP) or IFIT1/2 highlighted in ^23,25,28 are missing from both networks. Regarding the networks’ calibration, it allows for at most one computational state between measured time points. As the gaps between time points become larger (between 72 hours and 1 week for example), the freedom to incorporate more intermediate updates might be beneficial. However, they should be implemented carefully, considering the increased computation runtimes they imply. Furthermore, the nature of the data themselves and what they represent are of paramount importance in the Boolean network construction process. The study we used to provide the calibration data for the MVA-induced response network also contains cytokine concentration data, as well as antibody and T cell responses data ⁴⁷, which we could use in the future to enable our modeling approach to predict and interpret immune mechanisms more explicitly. We also limited our modeling efforts to gene expression and cell abundances in the blood. To broaden its scope and provide a more complex overview of the vaccine-induced responses, modeling the same parameters in lymphoid organs and other tissues, especially in the vicinity of the vaccine administration site, would be beneficial. These compartments could be interconnected within a single Boolean network or an agent-based modeling framework, as respectively illustrated for multi-cellular models in ¹²¹ and ¹²².

In order to validate such developments, increasingly rigorous evaluation efforts should be undertaken. A key challenge lies in the nature of the validation data itself. In therapeutics design, specifically for applications in oncology, increasing resources regarding “perturbomics”, i.e. the systematic perturbation of biological systems, provide a way to map cellular responses to stimuli and help refine predictive therapeutic models ^123,124. These perturbations, often chemical or genetic, are relatively straightforward to induce and measure, especially in vitro. In contrast, vaccinology presents a more complex landscape with immune responses that are local and systemic, which makes them difficult to replicate accurately outside of living organisms. Here, we leveraged the fact that mutant poxviruses have been extensively studied as vaccine vectors ²², often in a preclinical setting. A virtual evaluation of the use of pharmacological agents for host-targeted therapies, ideally stemming from drug repositioning ¹²⁵ together with MVA, as a strategy to enhance its immunogenicity or to model immune responses in groups of treated patients to personalize vaccines, could also be conducted. Ideally, simulated and experimentally perturbed responses should be integrated within a single study to strengthen the legitimacy of the modeling approach. This could be seen as the start of an iterative circle of model refinement and empirical validation, where each cycle enhances both the predictive power and biological relevance of the system. In this framework, simulations help generate hypotheses about immune dynamics, which can then be tested through targeted perturbations. The resulting data not only confirm or refute the model’s predictions but also feed back into the simulation pipeline, allowing for a recalibration of the model. Over time, this virtuous cyclic approach should give a deeper mechanistic understanding of vaccine-induced immunity, bridging the gap between in silica simulation and in vivo complexity. Moreover, this iterative cycle lays the groundwork for the emergence of digital twin technologies in vaccinology. Originally rooted in engineering, a digital twin refers to a virtual replica of a physical object or system, continuously updated with real-time data to mirror its behavior. Therefore, we argue that a sufficiently detailed and up-to-date computational model can serve as a digital twin of a biological process. Boolean modeling, notably, has proven particularly effective for precise and even individualized simulations ^65,126,127. Whereas emerging applications based on such a framework include personalized drug development and treatment optimization ^128–13°, we believe similar advances could be brought to the fields of vaccine design and optimization.

Supplementary figures

Supplementary Figure S1:

Supplementary Figure S2:

Supplementary Figure S3:

Supplementary Figure S4:

Supplementary Figure S5:

Supplementary Figure S6:

Supplementary Figure S7:

Data availability

Data complies with FAIR principles. Our findings could not have been be achieved without a computational tool we previously developed and published to model complex biological systems (doi: 10.1109/TCBB.2024.3456302), and we made the boolean models available for simulation, perturbation, and analysis on the CellCollective platform, a public deposition and execution environment.

Acknowledgements

We thank Nathalie Mantel (Vaccine Research and Development, Sanofi, Marcy 1’Etoile, France); Pierre Rosenbaum, Nicolas Tchitchek, Frédéric Martinon, and Roger Le Grand (IDMIT), Maria Chanzy (Dassault Systèmes) and Lina Cheraitia (IDMIT and Dassault Systèmes) for the valuable discussions.

Additional information

Funding

This work was supported by the Region Ile-De-France through the 2021 Paris Region PhD Program and Dassault Systèmes (for VD). The IDMIT infrastructure has been supported by the Programme Investissements d’Avenir (PIA), managed by the ANR under reference ANR-11-INBS-0008. The YF17D NHP study has received support from the Innovative Medicines Initiative Joint Undertaking (Biomarkers For Enhanced Vaccine Safety (BioVacSafe), IMI JU Grant Agreement no. 115308), resources of which are composed of financial contributions from the European Union’s Seventh Framework Programme (FP7/2007-2013) and EFPIA companies in kind contribution.

Author contributions

Original draft writing: VD; Conceptualization: VD, PC, MC, LN, ASB; Methodology: VD, PC, MC, LN, ASB; Review, editing, and approval: VD, PC, JGA, ME, MC, LN, ASB; Supervision: PC, MC, LN, ASB; Funding acquisition: PC, MC, LN, ASB.

Funding

Region Ile de France (2021 Paris Region PhD Program)

• Marine Ciantar

• Anne-Sophie Beignon

• Laurent Naudin

• Philippe Castera

• Vincent Deman